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Summary 

Simple  fluids  such  as  gasses  and  liquids  are  the  result  of  collisions  between  molecules. 
More  complex  fluids,  such  as  granular  flows  and  colloidal  suspensions  (non-Newtonian 
fluids),  result  from  the  more  complex  collision  (or  interaction)  behaviors  of  their 
constituent  particles.  In  this  project  we  have  demonstrated  that  collision  rules  can  be 
constructed  for  large  chunks  of  fluid  material  (eddies)  such  that  the  resulting  collective 
system  behaves  like  the  mean  (RANS)  flow  of  a  turbulent  fluid.  These  collision  rules 
are,  in  essence,  a  turbulence  model. 

The  project  has  demonstrated  its  three  primary  objectives.  First,  it  has  shown  that 
modeling  turbulent  flow  as  a  collection  of  colliding  (interacting)  objects  (eddies)  is  a 
theoretically  viable  approach.  Second,  the  project  has  shown  that  modeling  turbulence  in 
this  way  can  be  made  computationally  efficient  and  comparable  to  classic  Reynolds  stress 
transport  (RST)  models.  Finally,  the  collisional  approach  to  turbulence  modeling  has 
lead  to  some  insights  into  turbulence  and  turbulence  modeling  that  would  probably  not 
have  been  achieved  via  the  traditional  RST  approach. 

The  prediction  (or  modeling)  of  turbulent  fluid  flow  is  arguably  the  greatest  bottleneck  in 
the  Navy’s  ability  to  rapidly  design  innovative  devices  and  respond  to  environmental 
threats73.  While  this  research  does  not  address  a  specific  Navy  operational  issue,  it  has  an 
extremely  broad  impact  on  Navy  operations  and  the  Navy’s  ability  to  successfully 
execute  its  mission. 

Motivation:  Why  a  collision  model. 

The  collision  model  is  inspired  by  the  strong  analogy  between  granular  flows  and 
turbulent  fluid  flow.  An  example  of  a  granular  flow  might  be  Cheerios  in  a  factory  being 
piped  and  then  poured  into  boxes.  Similarly,  turbulent  flow  can  be  thought  of  as  a 
collection  of  eddies  that  interact  with  their  local  neighbors  as  they  are  piped  and  poured. 
In  both  cases,  the  particles  of  interest  are  of  roughly  the  same  size  as  the  pipes  and  mean 
flow  length  scales.  This  causes  the  resulting  flows  to  be  non-Newtonian.  The  behavior 
of  the  two  flows  is  not  identical  because  turbulent  eddies  have  a  range  of  length  scales 


and  Cheerios  (or  sand  grains,  or  many  other  granular  flows)  have  a  uniform  size.  In 
addition,  eddies  and  Cheerios  have  different  interaction  behaviors. 

The  traditional  approach  to  modeling  turbulence  or  non-Newtonian  fluids  is  to 
hypothesize  equations  for  unknown  stress  tensors  (in  turbulence  this  is  the  Reynolds 
stress  tensor).  Because,  the  particles  making  up  the  flow  are  roughly  the  same  size  as  the 
gradients  in  the  mean  flow  these  particles  respond  on  similar  timescales  as  the  mean 
flow.  This  means  that  algebraic  models  are  rarely  predictive,  and  evolution  equations  for 
the  stress  tensor  must  be  hypothesized.  In  turbulence,  these  are  the  Reynolds  stress 
transport  (RST)  equations.  Simpler  turbulence  models,  such  as  k-e,  are  simplifications  of 
the  RST  equations.  Similar  transport  equations  are  hypothesized  for  non-Newtonian 
flows,  and  many  turbulence  modeling  concepts  (realizability,  material  frame  indifference, 
tensor  consistency)  were  borrowed  directly  from  the  non-Newtonian  literature  at  this 
transport  equation  level. 

However,  it  has  long  been  recognized  in  the  non-Newtonian  fluid  community  that 
transport  equation  models  have  serious  limitations.  An  alternative  approach  is  to  model 
the  fluid  at  the  particle  collision  level  rather  than  using  a  transport  equation  for  the  stress. 
This  approach  is  more  versatile,  and  in  many  ways,  more  fundamental.  For  example, 
modeling  a  gas  as  particles  with  binary  elastic  hard  sphere  collisions  gives  the  Navier- 
Stokes  equations  and  the  perfect  gas  law  when  the  density  is  high,  as  well  as  the  correct 
gas  behavior  when  the  density  is  low  (where  Navier-Stokes  is  not  valid). 

Numerical  solution  of  collision  models. 


Once  a  certain  collision  behavior  has  been  hypothesized  there  are  three  very  different 
ways  to  solve  the  particle  system  numerically  and  obtain  a  prediction  of  the  fluid 
behavior.  The  most  straightforward  technique  is  the  ‘molecular  dynamics’  approach 
where  one  numerically  tracks  all  the  particles  in  the  domain,  and  performs  collisions 
when  they  occur.  This  approach  has  a  computational  cost  equivalent  to  large  eddy 
simulation  (LES)  and  is  not  considered  further.  The  other  two  approaches  note  that  one 
does  not  really  care  what  happens  to  individual  particles  but  only  what  happens  to 
particles  on  average.  The  quantity  of  interest  then  becomes  the  probability  density 
function  that  describes  the  probability  that  a  particle  (at  a  certain  place  and  time)  has  a 
certain  velocity.  The  evolution  of  the  probability  distribution  function,/,  obeys  the  exact 
equation 


df  df  df  df 

—+v— — ha 

dt  ‘  dx,  ‘  8v.  dt 


collisions 


(1) 


where  a,  is  the  acceleration  due  to  external  forces  (like  gravity)  and  the  right-hand  side 
describes  the  average  affect  of  the  collisions  on  the  PDF.  It  is  this  average  collision 
behavior  that  we  now  wish  the  models  to  predict. 

There  are  two  different  ways  to  solve  this  PDF  equation.  One  way  is  to  assume  the 
collision  model  has  a  Fokker-Planck  form  (see  equations  2  through  4).  Then  using  the 
equivalence  between  the  Fokker-Planck  equation  and  the  Langevin  equation  (Brownian 
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motion),  it  is  possible  to  construct  a  Lagrangian  particle  method.  This  is  the  approach 
extensively  researched  by  Pope  and  coworkers4’14,18’55.  The  Lagrangian  particles  move 
like  Brownian  dust  particles.  They  move  with  the  flow  and  are  randomly  perturbed  using 
a  prescription  given  by  the  model.  In  this  way  each  particle  is  independent  from  all  the 
others,  and  simply  interacts  with  the  average  of  all  the  other  particles.  This  is  less 
expensive  than  tracking  and  implementing  individual  collisions  (‘molecular  dynamics’ 
approach)  but  is  still  expensive  because  a  large  statistical  sample  of  particles  is  required. 

The  numerical  approach  used  in  this  project  was  to  solve  the  PDF  equation  using  a 
standard  Eularian  mesh  in  physical  space,  x,  as  well  as  in  velocity  space,  v.  Normally, 
this  approach  would  be  rejected  outright  since  10  mesh  points  in  each  direction  then 
requires  a  million  mesh  points  in  total  and  is  too  expensive.  The  resolution  to  this 
problem  is  to  use  an  extremely  coarse  mesh  in  the  velocity  space  (3  points  in  each 
direction).  This  means  we  are  solving  27  equations  for  each  point  in  space.  For 
comparison,  a  RST  model  solves  3  velocity,  1  pressure,  and  6  stress  equations  (10 
equations)  per  point  in  space.  However,  since  the  RST  equations  are  highly  coupled  and 
nonlinear,  and  the  PDF  equations  are  not,  the  solution  times  are  very  comparable. 

A  very  coarse  mesh  in  velocity  space  is  an  idea  borrowed  from  Lattice-Boltzmann 
methods  for  solving  the  Navier-Stokes  equations.  These  methods  solve  a  PDF  equation 
with  a  very  simple  collision  term  that  is  intended  to  give  Navier-Stokes  (Newtonian)  fluid 
behavior.  The  difference  here  is  that  we  solve  a  PDF  equation  with  a  much  more 
complex  collision  term,  which  results  in  RANS  behavior  for  the  fluid.  The  coarse  mesh 
is  acceptable  in  both  cases  because  the  interest  is  not  in  the  PDF  itself  but  in  its  lowest 
order  moments  -  the  mean  flow  and  the 
stresses.  These  low  order  moments 
can  be  reasonably  extracted  from  a  very 
coarse  approximation  of  the  PDF. 

Note  that  the  Langevin  approach  is 
equivalent  to  approximating  the  PDF 
with  a  random  sample,  and  a  large 
sample  is  needed  even  to  approximate 
the  low  order  moments  reasonably  well. 

The  Langevin  approach  is  slower 
because  it  provides  more  information 
(about  the  higher  order  moments). 

Unfortunately,  we  have  little  interest,  in 
engineering  turbulence  models,  in  the 
extra  information  the  Langevin  solution 
method  provides. 

Comparison  to  Lattice  Boltzmann  Methods. 

While  the  approach  taken  in  this  work  is  inspired  by  the  success  of  lattice-Boltzmann 
numerical  methods,  the  approach  is  significantly  different.  This  is  because  the  PDF 
governing  molecular  interactions  (Lattice-Boltmann)  has  a  variance  (width)  that  is  much 
larger  than  the  mean  and  which  is  essentially  constant  (related  to  the  speed  of  sound).  In 


Collision  Model 


PDF  Methods  Particle  Tracking 

( ‘molecular  dynamics  ’) 

▼  ^ 

Coarse  Discretization  Langevin  Equation 

(‘ lattice  methods’) 

Figure  1  Taxonomy  of  collision  model 
approaches. 
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contrast,  the  PDF  for  turbulence  has  a  variance  which  is  much  smaller  than  the  mean 
(turbulence  intensities  are  measured  in  percent),  and  which  can  vary  significantly  (in  time 
or  space).  This  is  illustrated  in  Figure  2. 


Figure  2.  Left:  typical  PDF  for  molecules.  Right:  typical  PDF  for  turbulence. 


To  capture  the  turbulence  PDF  with  only  three  points  it  was  necessary  to  have  a  moving 
adaptive  mesh  in  velocity  space.  In  order  to  avoid  losses  due  to  interpolating  one  mesh 
to  another  as  the  mesh  moves,  we  implemented  a  fully  conservative  scheme  in  which  the 
mesh  moves  continuously  in  time  (during  the  timestep).  This  uses  technology  previously 
developed  by  the  PI  for  moving  meshes  in  physical  space71.  In  actual  practice  the  PDF 
is  three-dimensional.  An  isosurface  of  an  actual  PDF  (the  50%  value)  is  shown  in  Figure 
3.  This  PDF  is  modeling  the  behavior  of  the  Le  Penven  et  al.  retum-to-isotropy  Case  IH 
>  0  experiment.  Note  the  fairly  large  changes  in  the  shape  and  size  of  the  distribution 
even  for  this  simple  experiment.  It  can  also  be  seen  in  this  figure  that  a  spherical  PDF 
corresponds  to  isotropic  turbulence. 


Figure  3.  Evolution  of  the  50%  isosurface  of  the  PDF  for  the  retum-to-isotropy 
experiment  of  Le  Penven  et  al.  (case  III  >  0). 


Theoretical  Analysis. 

Lundgren  (1967)  first  derived  the  exact  expression  for  the  collision  term  in  the  PDF 
evolution  equation  for  turbulence.  As  might  be  expected,  this  collision  term  can  not  be 
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expressed  solely  in  terms  of  the  PDF,  and  practical  solution  of  the  PDF  evolution 
equation  requires  modeling  the  collision  term.  In  this  work  we  have  focused  on 
generalizations  of  the  Fokker-Plank  collision  term.  In  its  simplest  form  this  collision 
term  has  the  form, 


dt 


collision 


=  3[a(v, 


(2) 


where  m(.  =  yjdy  is  the  mean  velocity  and  a  and  b  are  model  constants.  For  turbulence 
this  needs  to  be  generalized.  Pope  and  coworkers  use  the  form, 


£ 

dt 


collision 


& 


dx. 


(3) 


where  v';.  =  v;  -u;  is  the  fluctuating  velocity  and  the  first  term  (the  drift  term)  now  has  a 

matrix  model  parameter  Gtj ,  and  a  viscous  term  has  been  added  for  near  wall  (low  Re 

number)  calculations.  The  conversion  of  these  Fokker-Planck  models  to  a  Langevin 
equation  for  numerical  solution  dictates  that  the  diffusion  term  (with  b )  be  isotropic  and 
not  have  a  tensor  coefficient. 


In  this  project  we  analyzed  the  following  even  more  generalized  Fokker  Plank  model. 


£ 

dt 


collision 


5V; 


H„ 


jr 
A 


,1  a 


U 


3V; 


+  — 


dxn  jj  dt 


mesh 


3V; 


(4) 


The  last  term  on  the  right  hand  side  accounts  (exactly)  for  the  mesh  motion.  The  first 
three  terms  involve  model  tensors.  Sometimes,  these  tensors  are  isotropic  and  governed 
by  a  single  parameter.  The  viscous  terms  account  for  low  Reynolds  number  effects  and 
strong  inhomogeneity.  They  do  not  involve  any  additional  parameters  and  were  derived 
via  analysis  and  the  condition  that  the  model  be  exact  as  it  approaches  a  wall  (in  the 
laminar  sub  layer). 


Taking  the  zeroth  moment  of  the  collision  term  (Eqn  4)  gives  zero,  so  the  zeroth  moment 
of  the  PDF  equation  is  the  mass  conservation  equation.  The  first  velocity  moment  of  the 
modeled  PDF  equation  gives  the  momentum  equation. 


a»„  |  d(ulUn+Rm) 
dt  dxt 


(5) 


This  implies  that  the  acceleration  is  given  by  an  =-pn  +  (/*«,•,„ ),,  •  The  viscous 
contribution  to  this  acceleration  is  necessary  only  if  the  viscosity  is  not  constant. 
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Taking  the  moment  of  the  modeled  PDF  equation  with  respect  to  v’n  v'm  gives, 


r  +u  r M)=(CAtCA)+(H„+HJ 

dt  ox,  ox 


'  ~„K 


(6) 


where  Tnmi  =  jv>'m  v'(.  fdv  and  k=±Ru.  The  tensors  Gy,  Htj,  and  Jy  describe  the 

model.  Complex  dissipation  and  pressure-strain  models  can  be  implemented  via  these 
tensors.  From  this  analysis  it  is  clear  that  any  Reynolds  stress  transport  model  can  also 
be  implemented  as  a  generalized  Fokker-Planck  collision  model. 

The  equation  for  the  total  resolved  kinetic  energy,  Er  =  J-jV,v./Jv  -jRu ,  is 


dis  9  “i  3 

-^  +— [u,Er  +Uk(Rlk -fJU* )]  =  -(/>«,),,  +“nJRjn  -mAUi,j  +«;,)  + 


3x 


L^ll  (?) 


The  resolved  kinetic  energy  correctly  looses  energy  as  a  result  of  large  scale  dissipation, 
and  via  turbulence  production.  It  is  completely  specified  and  does  not  depend  on  the 
model  coefficients.  The  details  of  these  mathematical  analyses  are  presented  in 
Appendix  A. 


Practical  Implementation 

When  implementing  the  Fokker-Planck  collision  model  (Eqn.  4)  on  a  coarse  mesh,  it  is 
attractive  to  make  the  change  of  variables  f  =  In (/) .  If /is  close  to  Gaussian  (which  is 

expected)  then  /  will  be  close  to  parabolic.  This  parabola  can  be  accurately  resolved 
and  interpolated  by  the  three  points  available  in  our  scheme.  The  evolution  equation  for 
/  as  descirbied  in  Appendix  B  is, 


H„$-\  ,  H,  9f  $ 


u  9v,||  iJ  9v;  9v , 


+ 


a/ 


1 


_  «  ^  A 

a/  a/ 


,  fa  9/1(3  v,:  a/|] 

+  //  ,"[9x/9xJj|{/r  +  ^9vI.|J 


U 

a/ 1  9/  9/ 

^|j  +  //9^- 


(8) 


A 

While  there  are  more  terms  to  compute  in  this  version,  the  equation  for  /  is  much  more 
accurate  to  solve  numerically.  In  addition,  low  order  methods  and  simple  difference 
stencils  suffice  because  /  is  expected  to  be  very  close  to  quadratic. 
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The  models  for  the  tensors  Gy ,  Hi} ,  and  Jy  require  a  time  scale  to  be  dimensionally 

correct.  For  this  reason  an  additional  transport  equation  for  the  timescale  must  be 
included  in  the  model.  We  have  used  the  standard  epsilon  transport  equation  for  this 
purpose  since  it  is  very  commonly  used  in  RST  models  as  well. 

Summary  of  the  Model 

As  shown  in  Appendix  A,  the  above  PDF  equation  results  in  the  following  turbulence 
equation. 


^B-+^-u,Rmn  +  t^-vVv'  +  Rni  ^  +  Rmi  ^  =  (GmiRin  +  GniRim ) 

dt  ‘  mn  dxt  m  n '  m  a*,.  m'  0*,  v  mj  Jn  n]  jm) 

+  (H  +H  )  +  —  (7.  u  +J.  u  .)-2ju— 

\  mn  ±J-nm)  ^  \  in  m,i  im  n,i }  r*  ^  ^  J 


(9) 


The  triple  correlation  term, — v^v'v' ,  is  modeled  by 

3bt, 


a  m 

— vT — — 

dx,  dx. 


The  collision  model  is  given  by, 

o,j  =  c;,ss+c;,ivs +i(c;, 
»s=t  ecA+jKc^s,, 

J„=-jKC'rA 


(10) 

(11) 

(12) 


The  equivalent  Reynolds  stress  transport  equation  would  be, 

_a_ 

dx. 


dR„  a  „  a  a/?  /  „  „  \ 

— —  H - U:Rmn - vT  — —  +  ( u  .R.  +  u  ,R,m )  = 

av  ^  T  dx  '  ,J  1  ,J  1 


dt  dx , 


(c;2smJ+c;2wmJ)RJn+(c;2snJ + c;2waJ)Rjm +(c;2 -i)-Rm„ -jCplRmn  (13) 
+  +  4£C*„„S„.  +  -^-v^-  -  2ju^K  ^  R~‘  ^ 


^ p\  mn  1  3  J 


'p2  mn  dx.,  dx,  ^dx\K  l 


where  Si}  = — (mi  ;.  +  uj  t) ,  Wy  =—(ui  J  -  Ujy)  +  £ijk£lk ,  and  Q*  is  the  rotation  rate  in  a 
2  2 

non-inertial  frame  of  reference,  and  P  =  -Rnmun  m  is  the  turbulent  production  rate.  The 
model  parameters  are  given  by 


1.8- 
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C*  =-0.2F2  +  .006- 
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v 
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v  +  v , 
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where  F  =  -f-det^.  Ik)  is  the  two-component  parameter. 

The  transport  model  for  the  epsilon  equation  is  fairly  standard  and  is  given  by, 

T" +  u‘  =  T  (Ce'P  ~  Cel£  ^ +  T"  +  CeivT  )|—  ( 14) 

dt  ox.  K  dxi  dxi 

where 


CEi  =1.43 


Ce2  =11/6 


C£,  =0.33+0.5- 

£3  | 


Summary  of  Results 

The  PDF  collision  model  (using  a  coarse  moving  mesh)  was  tested  on  isotropic  decaying 
turbulence  at  different  Reynolds  numbers  and  rotation  rates.  This  flow  is  affected  by  the 
epsilon  equation  model,  but  not  by  the  PDF  model.  However,  it  is  a  good  test  of  the  PDF 
numerics,  since  the  exact  solutions  (for  a  given  epsilon  equation  model)  are  known.  It 
was  found  that  spacing  the  mesh  points  at  1.61  times  the  PDF  variance  was  on  optimal 
placement  that  almost  eliminates  numerical  errors  due  to  the  coarse  mesh  interpolation. 
The  moving  mesh  algorithm  moves  the  mesh  to  keep  it  at  this  position  even  as  the 
variance  changes  in  time. 

The  model  was  then  tested  on  anisotropic  decaying  turbulence.  Two  different 
experiments  (Le  Penven  et  al.,  and  Choi  &  Lumley)  and  five  different  data  sets  were  used 
to  evaluate  the  performance  of  the  model.  This  is  essentially  a  test  of  the  models  ability 
to  correctly  predict  slow  pressure-strain  or  retum-to-isotropy.  Some  very  interesting 
results  were  produced  by  this  study.  Whereas  all  RST  models  use  at  least  one  model 
constant  (and  frequently  two  or  more)  to  model  retum-to-isotropy,  the  collisional 
approach  to  this  problem  resulted  in  a  return  model  that  performs  well  and  has  no  model 
constants.  This  is  described  in  detail  in  Appendix  C,  which  has  been  submitted  for 
publication. 

Next  the  model  was  tested  in  a  variety  of  homogeneous  shear  flows.  They  key  to 
predicting  these  flows  correctly  is  in  the  modeling  of  the  fast  pressure-strain.  A  detailed 
study  of  the  fast  pressure-stain  model  parameters  was  performed.  These  parameters  are 
often  specified  as  constants,  however  they  can  be  functions  of  the  Reynolds  stress  and 
mean  flow  gradient  tensor  invariants.  A  wide  variety  of  DNS  and  experimental  cases 
were  used  to  back  out  the  exact  values  for  the  fast  pressure-strain  model  parameters  in 
these  experiments.  Both  rapidly  strained  and  slowly  strained  homogeneous  flows  were 
analyzed  in  this  way  to  obtain  the  parameter  values  in  both  limits.  With  this  data, 
models  for  the  fast  pressure-strain  parameters  were  proposed.  This  analysis  is  presented 
in  Appendix  D. 
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Finally,  the  model  was  implemented  and  tested  in  fully  developed  channel  flow.  The 
issue  here  is  to  correctly  account  for  inhomogeneity  and  low  Reynolds  number  effects. 
In  this  situation,  the  modeling  of  the  dissipation  tensor  requires  close  attention.  This  term 
dominates  near  the  wall  and  balances  viscous  diffusion.  A  model  which  accounts  for 
dissipation  in  inhomogeneous  flows  was  developed.  It  has  been  submitted  for 
publication,  and  is  presented  in  Appendix  E.  The  method  of  determining  the  model 
constants  is  shown  in  Appendix  F,  and  the  results  are  in  Appendices  G  and  H.  This  model 
for  the  dissipation  tensor  is  exact  in  regions  of  strong  inhomogeneity  and  involves  no 
model  parameters.  The  second  to  last  term  in  Eqn.  4  is  due  to  this  model.  The  fact  that 
the  model  is  exact  in  this  limit  is  important.  It  means  that  the  diffusion  is  exactly 
balanced  at  the  wall,  and  therefore  that  the  Reynolds  stresses  always  have  the  correct 
asymptotic  limits  near  a  wall.  This  means  that  elliptic  relaxation  approaches  are  not 
required.  In  addition,  computational  stability  is  significantly  enhanced  since  this  is  the 
region  where  Reynolds  stresses  are  close  to  becoming  unrealizable  due  to  the  numerics. 


Final  Conclusions  &  Future  Work 

This  project  has  allowed  us  to  demonstrate  that  collisional  models  are  a  viable  alternative 
to  RST  models.  In  one  instance,  we  have  even  been  able  to  remove  a  model  parameter 
due  to  insights  gained  from  this  viewpoint.  However,  it  is  also  clear  that  this  approach, 
as  it  stands,  has  most  of  the  same  difficulties  and  limitations  as  RST  models.  In 
particular: 

•  The  fast  pressure-strain  model  largely  dictates  the  model’s  performance  in  flows 
with  mean  flow  gradients  (most  flows).  There  are  a  large  number  of  constants  (a 
minimum  of  three),  practical  difficulty  in  modeling  equilibrium  and  rapid  limits 
within  a  single  model,  and  even  fundamental  problems  with  assuming  that  this 
term  can  be  modeled  in  terms  of  the  available  unknowns  (we  know  in  many  cases 
that  it  can  not).  These  difficulties  are  common  to  both  RST  models  and  the 
current  collision  model. 

•  The  scale  (or  epsilon)  transport  equation  is  the  other  source  of  significant  error 
and  parameterization  (many  constants).  In  the  simple  form  we  have  so  far 
studied,  the  collision  model  in  no  way  address  this  issue.  The  same  scale 
transport  equation  is  used  as  in  RST  models. 

•  Finally,  although  we  have  used  Lattice-Boltzmann  discretization  ideas,  the 
implementation  of  these  models  is  not  as  computationally  efficient  as  an  analysis 
of  Lattice-Boltzmann  methods  might  first  suggest.  The  fact  that  a  moving 
adaptive  mesh  is  required  for  collisional  turbulence  simulations  means  that  most 
of  the  speed  enhancing  tricks  of  Lattice-Boltzmann  methods  can  no  longer  be 
applied.  The  method  is  therefore  computationally  also  comparable  to  RST 
models. 

One  important  distinction  between  RST  models  and  the  collisional  modeling  approach  is 
that  many  of  these  deficiencies  can  actually  be  fixed  within  the  collisional  framework. 
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The  key  is  to  include  more  information  in  the  PDF  than  just  the  velocity.  By  assuming 
that  eddies  have  shape  (or  more  specifically,  orientation)  as  well  as  velocity,  two  major 
problems  can  be  addressed.  First,  the  fast  pressure-strain  term  can  now  be  implemented 
exactly.  Since  the  production  term  is  already  exact,  this  means  that  the  influence  of  the 
mean  flow  on  the  turbulence  is  now  captured  exactly  (up  to  numerical  implementation 
errors),  only  nonlinear  (turbulence-turbulence)  interactions  must  be  modeled.  Secondly, 
eddy  shape  provides  a  length  scale  measure,  so  a  separate  scale  transport  equation  is  no 
longer  necessary.  Including  orientation  into  the  collisonal  model  will  increase  the  cost 
by  an  order  of  magnitude.  However,  the  number  of  model  parameters  will  also  be  vastly 
reduced.  If  the  resulting  model  is  highly  predictive  then  the  additional  cost  may  be 
warranted  and  will  still  be  far  below  the  cost  of  the  next  alternative  (large  eddy 
simulation). 
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Appendix  A:  Moments  of  the  PDF  Equation 


In  theory,  the  PDF  contains  enough  information  to  calculate  all  single  point  statistics 
involving  the  velocity.  This  means  that  from  a  particular  collision  model  we  can  derive 
the  resulting  mass,  momentum,  total  kinetic  energy,  and  Reynolds  stress  equations.  In 
this  appendix  a  very  general  model  form  is  assumed,  and  the  resulting  equations 
determined. 

The  general  PDF  evolution  equation  is  given  by. 
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(J V  +Wj) 
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J  dv, 
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■)Xj 

Hrt[K)  .1 
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Momentum: 
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1  dt  dtJ  dt 

(A.  3) 

Turbulence: 

J 

[v'v'a/  dvJR™ 

1  m "  dt  dt 

(A.4) 

Energy: 
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f  ^ 

(A.5) 

2. 


_  _  r  df  d  r  dut 

Mass:  Yv—dv  =  —  \vjdv  =  — L 

J  oxi  oxi  J  oxi 

Momentum:  \vivk^7~dv  =— \vivkfdv 

J  axi  dxi  J 

Writing  v  in  terms  of  v  gives: 


=  ^\(v'  +  ui)(v'k+uk)fdv 
=  j-  {(vX  +  v'uk  +  u/k  +  uku.)fdv 


(A.6) 
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Turbulence: 


dRik  ,  du,uk 
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Since  J«,v^/dv  =  0  we  get  the  following. 
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Convection  Production 

Turbulent  Transport 


Energy:  j  fw^dv  =  \  JviVj 


9x 


_«  ^ 
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=y[l(M,+v0^v*]=y-  iu,vkvk  +lv'(w,+v;)2|j 


Since  Ji^v'/dv  =  0  we  get  the  following. 

=  Ar„  777  j.  1.77777' 

9x,  I 
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9x 
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Mass: 


(A.22) 


\jr-[{ai  +  Guv'j)f]dV=  1  (ai  +  GijVj)fhids^O 

ovi 

By  using  Gauss’s  divergence  theorem,  we  can  convert  this  volume  integral  into  a 
surface  integral,  where  n  is  the  normal  vector  to  the  cell  face.  Evaluation  this 
integral  over  all  velocity  space  means  this  term  is  zero,  since  the  probability 
density  function  goes  to  zero  as  v  foes  to  infinity.  . 

Momentum:  jv*  +  Gyv'j )  /  J  dv 

OVi 

=  (“- + G/,)-f]‘K  -  Jk  <A-23) 

=  j  v,  (a,  +  GyVj )  fiifdS  -  at  Jfdv  -  Gy  jv'jfdv  (A.24) 

9Q— 


The  first  term  is  zero  by  Gauss’s  Theorem  and  the  last  because  the  average  of  the 
fluctuations  is  zero.  jv'mfdv  =  0 . 


=  -at 


Turbulence: 


Kv^[(«,+^v;)/]^ 


Writing  v  in  terms  of  v  gives: 
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By  Gauss’s  Theorem  we  can  eliminate  the  first  term. 
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By  Gauss’s  Theorem  we  can  eliminate  the  first  and  third  terms. 
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Since  jukv'jfdv  =  0  we  get  the  following. 
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By  Gauss’s  Theorem  we  can  eliminate  the  first  term. 

=-Kv,fU 

Again,  from  Gauss,  we  see  the  first  term  is  zero. 
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Turbulence: 


r  3  df  3  3  r  - _ 

Jdxi  dxi  axi  dxi  J 

r  3  3/  3  3  r  3  duk 

Jv*  =  3“^T~  Jv*^v  =  3“^^ 

J  dxt.  3x(  dXj  dxt  J  dxi  dx{ 

t ,  ,  3  3  f  j- 

lv  v — il-^—dv 
axj  dXj 

3/.l3v'v' 


r  3 

'  a/  ,1 

J3x 

w/—  - 

L  a*j  |_ 

(A.46) 

(A.47) 


(A.48) 


We  note  that  — ^  =  -v' w„  .  -  v'«  , 

w  /j, /  n  m,i  7 

dXj 


=  1: 


3x 


/  /  a/,1 

V  V  ll—\  + 

m  tu  1 1 

3xJ 


3 f  .1  /  /  /  \j- 

.  d*.j 


r  3  I"  3 v'  v  f  3v' v  1  3/  /  ,  /  \ 

+v*“-')dv 


(A.49) 

(A.50) 


=  /v' v'/3v  +  j-  j///(v>ni  +  v>ffl .  )dv  +  //  )Jv 


3x,  3x, 


3x, 


3x 


The  middle  term  is  zero  due  to  the  fact  that  jv'mfdv  =  0 
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Again  the  middle  term  is  zero. 
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Using  Gauss’s  theorem, 
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Turbulence:  {4^(7..  +fiuJJ)^J-dv  =  {J..+ ^Uji)\vmvnj-^-dv  (A.65) 
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Gauss’s  divergence  theorem  eliminates  the  first  term. 
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Second  term  is  zero 

3  rM 


=-^'arf 

=-2juKt 

Energy:  \\vkvk{nKt 


dxli  K 

O 

K  ), 

d_d_ 

o be,  dvt 


+^-?-dv 
K 


XT 

*!b 


Vli 


d_d_ 

dVj  dxl 


vkvk 


dx, 


-/L 


AH 


K 


By  Gauss’s  Theorem  we  can  eliminate  the  first  term 


Ml 


£ 


If 


Since  jv'/Jv  =0  we  get 


Vf 


=  0 

So,  in  summary. 


Mass: 


0  =  - 


9m, 

dx, 


n .  .  du.  dR..  duu. 

Momentum:  +  -H*-  H — -  k 


9  du. 


~“k=^M- 


dt  dx.  dxt  *  dx; r  dXj 

Turbulence: 

9i?„„  d  _  d  „  du  du  /  \ 

+^r«,«- + 37  w, + K  af + R.,  = (Gw«* + ) 


(A.  8  6) 

(A.  87) 
(A.  8  8) 

(A.89) 

(A.90) 

(A.91) 

(A.92) 

(A.93) 

(A.94) 

(A.95) 
(A.  96) 


dt  dx. 


dXj 

+  (Hmn+Hnm)+— ju^^-(j.u  +J.u  )-2u—^—^ 

\  mn  nm )  ^  r*  m  m,t  im  n,i )  r*' 


V  K  j, 


(A.97) 


Energy: 

f)/?  3  r  _  ^ 

a7 +V’Lu'jB + UkRik  +  *  w*  J = fl‘M< + GA +/4  +  57/*— ) 


9m, 


9x 
(A.98) 


8 


Appendix  B :  NinnoMSoliitionofthe  PDFEquatkxi 


PDF  Equation: 
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‘F  hat’  Equation: 


By  defining  ‘F  hat’  as  the  natural  log  of  the  PDF  function,  we  arrive  at  an  equivalent 
evolution  equation,  which  looks  like  the  following 
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The  resulting  evolution  equation  for  ‘F  hat’  is: 
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Appendix  C:  PDF  Models  for  Return- to-Isotropy 


1.  Introduction 

In  most  turbulent  flows  of  interest  the  turbulent  velocity  fluctuations  are  anisotropic,  that 
is  they  differ  in  magnitude  depending  on  their  orientation.  One  aspect  of  Reynolds  stress 
transport  models  (and  other  more  advanced  models)  that  distinguishes  them  from  simple 
two-equation  transport  models  like  k-e  is  their  ability  to  more  accurately  model 
turbulence  anisotropy.  The  degree  of  anisotropy  is  important  because  it  can  directly 
impact  how  turbulence  affects  the  mean  flow. 

In  the  absence  of  any  driving  mechanism  anisotropic  turbulent  flows  tend  to  return  to  an 
isotropic  state  (the  state  of  least  order).  This  nonlinear  process  is  often  called  retum-to- 
isotropy.  It  was  identified  early  on  in  the  development  of  Reynolds  stress  transport 
models  and  first  modeled  by  Rotta1.  Since  that  time,  the  retum-to-isotropy  process  has 
been  extensively  investigated  and  modeled2'10. 

The  retum-to-isotropy  problem  is  of  significant  theoretical  interest  in  the  theory  of 
turbulence  because  it  is  entirely  due  to  nonlinear  interactions.  The  return  process  is 
irreversible,  and  is  the  mathematical  consequence  of  the  averaging  process.  As  a  result 
we  know  that  the  return  process  must  be  modeled,  and  that  there  can  be  no  exact 
representation  for  this  process  in  terms  of  Reynolds  averaged  variables.  Existing  models 
for  retum-to-isotropy  tend  to  make  extensive  use  of  mathematical  concepts,  such  as  the 
Cayley-Hamilton  theorem,  realizability,  Taylor  series  expansions,  and  fixed  point 
analysis.  The  resulting  models  invariably  have  a  least  one  model  ‘constant’  that  must  be 
set  via  experiments. 

In  this  work,  we  are  interested  in  deriving  models  for  the  return  process  based  on  physical 
ideas  as  well  as  the  commonly  used  mathematical  tools.  We  make  the  assumption  that 
turbulence  behaves  as  a  kinetic  process,  and  that  kinetic  models  of  turbulence  may  lead  to 
some  useful  insights  about  the  return  process.  The  advantage  of  this  approach  is  that  the 
resulting  models  can  be  made  to  be  free  of  any  tunable  model  constants. 

In  Section  2,  the  classic  Reynolds  stress  transport  equation  approach  to  modeling  retum- 
to-isotropy  is  briefly  reviewed.  We  use  these  classic  results  as  a  reference  since  this  is 
the  approach  that  is  most  widely  understood  by  most  readers.  In  Section  3  we  consider 
retum-to-isotropy  from  the  perspective  of  the  BGK11  approximation  to  the  Boltzmann 
equation.  Classic  linear  return  models  result  from  this  kinetic  equation.  The  deficiencies 
of  the  BGK  approach  are  largely  solved  by  two  parameter-free  relaxation  collision 
models  developed  and  tested  in  Section  4.  Section  5  investigates  the  predictive 
performance  of  these  models  for  five  different  experimental  cases.  The  relaxation  model 
is  extended  in  Section  6  to  enable  any  desired  Reynolds  stress  return  behavior,  and 
another  parameter  free  model  is  proposed  that  has  some  unique  and  interesting  properties. 
Section  7  explores  the  implications  and  connections  to  the  Fokker-Planck  collision 
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model,  and  the  results  are  discussed  in  Section  8,  where  some  speculation  is  presented  as 
to  what  these  kinetic  models  imply  about  turbulent  eddy  interactions.  , 


2.  Reynolds  Stress  Transport  Models 


In  the  absence  of  any  mean  flow  the  evolution  of  the  Reynolds  stress  tensor,  ,  in 
homogeneous  but  anisotropic  turbulence  evolves  according  to  the  equation 


3f 


“2 VWj,k  +P(u’,j  +“}./> 


(C.l) 


The  first  term  on  the  right-hand  side  is  the  dissipation  rate  tensor  and  the  second  term  is 
the  slow  pressure-strain.  The  pressure-strain  is  considered  ‘slow’  in  this  situation 
because  the  pressure  in  this  term  depends  only  on  the  turbulence  not  on  the  ‘rapid’  mean 
flow  velocity  gradients  (since  there  are  none  in  this  situation).  Both  terms  require 
modeling.  However,  one  half  the  trace  of  the  dissipation  tensor  is  the  dissipation  rate, 

e  =  vu'i  ku\  k ,  which  is  assumed  to  be  known  (given  by  another  transport  equation  model), 
and  the  trace  of  the  pressure-strain  term  is  zero  in  incompressible  flows. 

The  most  common  modeling  approach  is  to  assume  that  the  dissipation  tensor  is  close  to 
isotropic.  If  small  anisotropy  in  the  dissipation  tensor  exists  then  it  is  then  included  with 
the  pressure-strain  model.  The  slow  pressure-strain  and  anisotropic  dissipation  are  then 
collectively  modeled  as  a  ‘retum-to-isotropy’  term.  There  are  reasons  to  suggest  that 
modeling  dissipation  anisotropy  and  slow  pressure-strain  separately  is  advantageous,12 ,13 
but  for  simplicity  we  retain  the  ‘collective’  approach  described  above.  The  simplest 
model  (due  to  Rotta)  is  that  the  retum-to-isotropy  term  is  proportional  to  the  Reynolds 
stress  anisotropy.  This  gives  a  Reynolds  stress  transport  model  of  the  form, 


dt 


(C.2) 


The  ‘retum-to-isotropy’  term  will  tend  to  drive  the  Reynolds  stress  tensor  towards  an 
isotropic  state  as  time  proceeds.  The  rate  at  which  this  happens  is  governed  by  the  Rotta 

A 

constant,  CR .  This  return  model  is  the  simplest  possible,  and  is  linear  in  the  Reynolds 
stress  anisotropy,  atj  =(-jf-f  <5jy)-  Equation  C.2  appears  to  imply  retum-to-isotropy  for 
any  positive  value  of  CR .  In  fact  this  is  not  the  case,  CR  must  be  greater  than  1.  To  see 
this  we  look  at  the  evolution  equation  for  which  should  tend  towards  •£  <5^. . 


<HVf)  %  3 K 

dt  K  dt  dt 


(C.3) 


The  isotropic  dissipation  actually  causes  the  Reynolds  stress  tensor  to  move  away  from 
isotropy  which  must  be  counteracted  by  the  return  term.  CR  is  actually  a  parameter,  not  a 
strict  constant,  that  can  be  (and  often  is)  a  function  of  the  Reynolds  stress  invariants  and 
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turbulent  Reynolds  number.  Due  to  the  strict  requirement  described  above  the  splitting 

/\ 

CR  =  1  +  CR  is  useful.  This  gives  a  model  equation  of  the  form, 

^.=-e^-CRe^S„)  (C.4) 

where  CR  >  0.  Typical  values  for  CR  lie  between  0.5  and  1.0  (Durbin)14.  Launder,  Reece 
and  Rodi15  suggest  a  value  of  0.8.  No  return  to  isotropy  is  the  case  of  CR=  0. 
Physically,  the  no  return  limit  appears  to  occur  at  low  Reynolds  numbers.  In  addition, 
the  no-retum  limit  is  often  enforced  in  the  two  component  limit  (where  one  of  the 
Reynolds  stress  diagonals  goes  to  zero  faster  than  the  others,  such  as  near  walls).  For 
this  reason  CR  is  often  not  a  constant  but  is  actually  a  parameter  that  depends  on  the 
turbulent  Reynolds  number  and  Reynolds  stress  invariants616. 

It  is  helpful  to  propose  a  general  model  for  the  Reynolds  stress  evolution, 

£ = -£(rL  Kj + n,m  Rmi )  (c.5) 

where  the  dimensionless  Fly  is  some  as  yet  unspecified  model.  Expanding  this  model  as 
Uij  =  Sy +tlu  gives 

(c.6) 

A 

so  it  is  clear  that  II, y  is  the  return  part  of  the  model.  The  trace  of  the  last  term  should  be 

A 

zero,  so  we  have  a  single  constraint  on  the  model,  II, y  Rp  =  0 .  It  is  not  necessary  that 

A. 

II, y  be  symmetric.  The  explicit  inclusion  of  the  Reynolds  stress  in  Eq.  (C.5)  means  that 

this  general  model  can  be  strongly  realizable  (Schumann17,  Lumley2)  if  n,y  is  finite.  If 
one  component  of  the  turbulence  goes  to  zero  then  this  model  will  also  make  the  time 
derivative  of  that  component  go  to  zero.  However,  in  the  unusual  circumstance  that  n(> 
becomes  singular  (goes  to  infinity)  this  model  can  potentially  violate  strong  realizability. 
The  classic  linear  return  model  described  above  is  given  by  fly  =  CR(Sy  -jKR^1) .  This 
model  becomes  singular  in  the  two  component  limit  (because  of  the  Reynolds  stress 
inverse).  The  classic  linear  model  satisfies  weak  realizability18  if  CR  >0,  but  for  the 

linear  model  to  satisfy  strong  realizability  CR  must  go  to  zero  in  the  two  component 
limit. 

Slightly  more  complex  nonlinear  models  for  retum-to-isotropy  have  the  general  form, 

■f" it =  Cr  ( ay )  +  CN  ( ajkakj  —  -j- j  (C.7) 
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Cubic  and  higher  order  nonlinear  models  can  also  be  represented  by  this  quadratically 
nonlinear  model  due  to  the  Cayley-Hamilton  theorem.  Sarkar  and  Speziale3  suggest 
values  of  CR  =  0.7  and  CN  =  1 .05 . 

The  realizability  conditions  are  clearer  when  this  model  is  written  in  terms  of  the 
Reynolds  stresses, 

i  =  (C.8) 

Pre  and  post-multiplying  this  expression  by  the  eigenvector  matrix  Q  diagonalizes  the 
Reynolds  stress  tensor  ( QrRQ  =  D ),  so 

tffQ— fD-{C,-C,[i^-f])f(D-in)+C„^(DD-^D)  (C.9) 

since  Qr^Q  =  |p+D(^Q)-(-^Q)D  the  off  diagonal  evolution  is  trivial,  and  the 
diagonal  components  individually  satisfy  the  right-hand  side  of  Eq.  (C.9).  Weak 
realizability  is  satisfied  as  long  as  CR  -  CN  [-^r-  -|]  >  0 .  Strong  realizability  requires 

equality  on  this  previous  expression  and  1  +  CN  [-^pr]  ^  0 .  The  quantity  -^r-  appears 
frequently  and  is  related  to  the  second  invariant  of  the  anisotropy  tensor  via 


The  model  expression  for  the  nonlinear  return  model  is 

n„  =  S, .  +  ( C,  -  C„  [*£•  - 1  ] }  (*,  -  1KR7' ) -  C„  (i  - -  6, )  (C.  1 0) 

The  singularity  due  to  K?  is  weakly  realizable  as  long  as  the  leading  coefficient  is 

positive.  It  is  strongly  realizable  if  this  leading  coefficient  is  zero  in  the  2-component 
limit  and  the  coefficient  of  Stj  is  positive. 


3.  BGK  Collision  Models 


In  homogeneous  turbulence  in  the  absence  of  mean  accelerations  or  mean  pressure 
gradients  the  evolution  equation  for  the  velocity  probability  density  function  (PDF)  is 


M.-4L 

collisions 


(C.11) 


This  equation  governs  the  decay  of  anisotropic  homogeneous  turbulence,  which  is  the 
focus  of  this  work.  One  of  the  simplest  collision  models  is  a  relaxation  of  the  PDF  to 
some  known  equilibrium  form. 

I •  =  -f'CMr(/-/")  (C.  12) 
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where  the  constant  CBGK(x,t )  might  be  a  function  of  position  and  time  but  is  not  a 
function  of  the  velocity.  This  model  is  similar  to  the  BGK  approximation  for  collisions 
used  in  Lattice-Boltzmann  methods.  The  constant  CBGK  should  always  be  greater  than  0 
for  a  well-posed  method.  Unlike  molecules,  turbulence  particles  do  not  conserve  kinetic 
energy  when  they  collide,  so  the  form  of  feq ,  the  equilibrium  target  distribution,  must  be 
slightly  different  from  classical  theory.  If  we  take  the  target  distribution  to  be 

feg(K)  =  (j^kyV2e~ht^  (C.13) 

where  0  <  K  <  K ,  then  (as  shown  in  Appendix  I)  mass  and  momentum  are  conserved  and 
turbulent  kinetic  energy  obeys  the  equation,  =  — j  CBGK  ( K  -  K) .  This  implies  that 

CBGK  =  ,  and  the  dissipating  collision  model  is 

(cm) 

This  is  a  model  in  which  the  PDF  relaxes 
towards  a  spherical  Gaussian  PDF  with 
less  turbulent  kinetic  energy  (see  Figure 
1).  Those  portions  of  the  PDF  which  are 
farthest  from  the  target  spherical 
distribution  decay  faster  than  those 
portions  of  the  PDF  which  are  closer  to 
the  target.  Figure  1;  BGK  relaxation 

model .  Solid  line  represents 

The  equivalent  Reynolds  stress  transport  an  isocontour  for  an 

.  .  ,  .  -  .  .  .  t  anisotropic  PDF.  Dashed 

equation  is  obtained  by  multiplying  by  line  is  the  spherical  target 

v'jV'j  and  integrating  over  all  velocities. 

This  is  shown  in  Appendix  I,  and  results  in  the  following  equation. 

(C.15) 

In  terms  of  II,  this  model  is  II, •,  =£^L(Si]-^KR~j  ) .  Which  is  identical  to  the  classic 
return  model  if  CR  =  ,  or  equivalently  K  =  K-^ y.  This  implies  the  relation 

A 

CBGK  =\  +  CR=CR  between  the  BGK  relaxation  constant  and  the  Rotta  constant. 

From  this  analysis  it  can  be  seen  that  there  is  no  return  to  isotropy  if  CR  =0  (or  K  =  0). 

A 

Under  the  condition  K  =  0,  feq  becomes  a  delta  function.  This  observation  suggests  an 
alternative  model  of  the  form, 

l  =  -f(/-^(v’))-fQ(/-/e9(K))  (C.16) 
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The  first  term  (involving  a  delta  function)  produces  pure  decay  and  the  second  produces 
return  to  isotropy  with  no  decay  (relaxation  to  a  spherical  PDF  of  the  same  energy).  This 
two-part  model  has  been  proposed  by  Degond  &  Lemou10. 

While  both  C.14  and  C.  16  result  in  an  identical  equation  for  the  Reynolds  stress  evolution 
(the  classic  linear  Rotta  model),  the  models  themselves  are  not  identical.  Differences 
exist  in  the  evolution  of  the  higher  turbulence  moments.  The  model  given  by  Eq.  (C.16) 
will  tend  to  produce  a  spike  in  the  PDF  around  its  mean  value  (due  to  the  delta  function). 
Eq.  (C.14)  has  a  smoother  influence  on  the  PDF  in  general  but  will  also  produce  a  spike 
if  CR  goes  to  zero  (in  the  two-component  or  low  Reynolds  number  limits). 

Neither  model  has  the  ellipsoidal  (Eq.  23)  or  spherical  (Eq.  18)  Gaussian  as  a  solution. 
This  implies  that  even  if  the  turbulence  starts  with  a  Gaussian  PDF  it  does  not  stay 
Gaussian.  It  is  not  a  strict  fact  that  turbulence  should  be  Gaussian.  Certainly  under  the 
influence  of  inhomogeneity  we  know  it  is  not  Gaussian  at  all.  Even  in  homogeneous 
turbulence  the  tails  of  the  PDF  are  not  expected  to  be  Gaussian.  However,  statistical 
arguments  based  on  the  central  limit  theorem  would  suggest  that  decaying  homogeneous 
turbulence  ought  to  be  close  to  Gaussian  or  at  least  evolve  in  that  direction  for  most  of 
the  core  portion  of  the  PDF.  Experiments  (Tavoularis  &  Corrsin)19  of  homogeneous 
turbulent  shear  flows  support  the  hypothesis  that  homogeneous  turbulence  (even  when 
sheared)  has  a  central  core  that  is  closely  approximated  by  an  elliptical  Gaussian  PDF 
(sometimes  called  a  trivariate  normal  distribution). 


4.  Relaxation  Collision  Models 

A  more  general  form  than  the  BGK  model  (Eq.  C.12)  for  collisions  is  the  linear 
relaxation  model, 

| =g(y)-h(y)f  (C.17) 

where  g(v)  >  Oand  h(\)  >  0  are  some  positive  functions  of  the  velocity  (and  possibly 
position  and  time  as  well).  The  positivity  requirements  keep  the  governing  equation 
stable  and  the  probability  always  greater  than  zero. 

In  addition,  the  model  should  conserve  the  total  probability  (or  mass),  so  that 
J  g(v)d\  =  |  h(\)fdv ,  and  it  should  not  cause  any  mean  flow  to  be  created,  implying 

Jv n[g  ~ Hf]d\  =  0 .  Finally  the  model  should  dissipate  energy  at  the  correct  rate, 
]^[hf-g]dv  =  £. 

One  way  to  determine  a  suitable  choice  for  the  model  functions  is  to  insert  a  desired 
solution  for  the  PDF  function  /  and  then  derive  the  parameters  from  Eq.  (C.17).  In 

isotropic  decaying  turbulence  there  is  evidence  that  the  core  of  the  PDF  is  very  close  to  a 
Gaussian  and  retains  this  shape  during  the  decay  process  (Yeung  and  Pope)  .  If  we 
assume  the  PDF  equation  (17)  has  a  Gaussian  solution, 
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(C.18) 


/(v,o=(j^r'v 


where  v  \  =  vn  -  un  and  un  is  the  mean  velocity,  then  taking  the  time  derivative  gives, 


a t 


=  (jxK)~V2e~h^(l 


V’nv'r  \  3  g 
2K  >  2  K 


(C.19) 


Comparing  with  Eq.  (C.17)  suggests  that  a  suitable  choice  for  the  model  functions  is 
g(v)  =  feq{v)jj  and  h(v)  •  Actually,  these  functions  do  not  conserve 

momentum  or  dissipate  energy  at  the  correct  rate.  They  must  be  generalized  slightly  to, 


*(v)  =  CM#( UK) 


- 

e  4K 


h(y)  =  C, 


u 


M  IK  l2K+(u-u)7] 


(C.20) 


where  we  expect  CM  ->  1 ,  K  ->  K ,  v’ ,  v’,.  -4  v’,.  when  the  PDF  approaches  a  spherical 

Gaussian  (Eq.  (C.18)).  Conservation  of  mass  is  already  satisfied.  Conservation  of 
momentum  implies  a  relationship  exists  between  the  hat  and  tilde  velocities, 

(«p  -  ■ up )  [2K  +  (u  -  u)2  ]  =  2 Rip  ( u.t  -  u, ) +  Jv  pv  \  v  \  fdv  (C.2 1 ) 

This  implies  that  either  up  or  up  can  be  specified  arbitrarily  and  then  the  other 
determined  by  Eq.  (C.21).  The  two  simplest  choices  are  up  =  up  which  implies 

«,  =ui+^t  Jv ‘p v v fdv ,  and  iip  =  up  which  implies  u,  =m,. Jv',v'nv fdv .  In 
either  case,  if  the  PDF  is  symmetric  then  the  odd  order  integral  is  zero  and  up=up=up. 
Since  by  definition  v,.  =  w,+v  =  w,+v  ’.  =  m(.  +  v  \ ,  this  also  implies  v =  v =  v as  well. 

Therefore  the  hat  and  tilde  quantities  in  Eq.  (C.20)  can  be  viewed  as  a  small  perturbation 
imposed  when  the  PDF  is  skewed  (not  symmetric),  and  are  largely  a  formal  technicality 
to  enforce  conservation  of  momentum. 


Conservation  of  energy  imposes  a  relation  between  CM  and  K/K  (Appendix  J). 

[f+^(u-u)2+^][2K  +  (u-u)2] 

=  2*  Jv>Vv ’/ V \  fdy  +  {-^iA  Jv 'p V v fdv +  (u- u)2 


(C.22) 


If/is  symmetric  this  simplifies  considerably  to  -f+^f \v'pv'pv\ v\fdv.  If/is 
an  elliptic  Gaussian  given  by 


/  =  I(2*r)  det(/?  )]”  e  * 


~if2  \Rnmv  nv  m 


(C.23) 


then  the  integral  can  be  evaluated  and  is  4K2  +2RnmRmn  (Appendix  K).  Then 
f  =  or  perhaps  even  more  informatively  \CU  -  (l+-^r- -•§-)"’.  The 
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relaxation  model  therefore  has  one  free  parameter  (either  CM  or  K/ K).  Both  of  these 
parameters  should  go  to  1  when  the  turbulence  is  isotropic  (i.e.  when  /  is  a  spherical 
Gaussian).  Since  -» |  in  isotropic  turbulence,  forcing  one  of  these  conditions  is 

sufficient  to  guarantee  the  other. 

The  derivation  of  the  equivalent  Reynolds  stress  equation  is  given  in  Appendix  J.  The 
result  is  that  Eq.  (C.20)  is  equivalent  to 


dRij  _ £ _ /  Rij  ,  RinRnj  2  K  S\  \ 

dt  “  LMat.il  \  *  K 2  3  K  uij  ) 

21 K2  *U 


(C.24) 


if  an  elliptic  Gaussian  is  assumed  for  the  PDF.  Eq.  (C.24)  in  turn  implies  the  return 
parameters 


A_ 

RmnRnm 

3 

2  K2 

31 

1+ 

Rmn^nm 

A1 

2  K* 

^lj 

or,  in  terms  of  Cm, 
C  -1C  -1 

2  M  1 


and  CN  r  jwji*  *i 

Ll+  2*2  *y 

and  CN-—-^CM 


(C.25) 


(C.26) 


Note  that  this  model,  and  the  other  models  derived  in  this  work,  tend  to  imply  that  Cn  is 
less  than  zero.  In  contrast,  the  widely  used  nonlinear  model  of  Sarkar  and  Speziale3  has  a 
positive  value  for  this  constant.  The  implications  of  this  difference  are  examined  in 
detail  in  Section  8. 


Various  choices  of  Cm  are  possible.  The  simple  choice  CM  =1  leads  to  CR=j,  and 
CN  =-f .  These  values  produce  a  model  which  is  very  similar  to  the  two  models 
examined  in  detail  below. 


The  equally  simple  choice  -§-  =  1  implies 


3  R, 


2K 


2K 


iRnm 


md  c»  =-«. 


(C.27) 


This  choice  of  4  implies  the  ‘target’  distribution  has  the  same  energy  as  the  PDF  but  a 
spherical  shape.  The  performance  of  this  model  is  shown  in  Section  5  and  is  referred  to 
as  Model- 1.  The  realizability  condition,  CR  -CJ^r~i]  =tjt~  >  indicates  that  Model- 
1  is  weakly  realizable. 

In  general  the  realizability  condition  for  these  relaxation  models  is 
CR-C +  ~|)  so  choices  where  J-  vanish  in  the  2-component 

limit  will  satisfy  the  strong  realizability  condition.  The  quantity  F  =  det(Rij)/(jK)3  is  1 
in  isotropic  turbulence  and  0  in  the  two  component  limit.  The  choice  jf  =  F  means  that 


8 


(C.28) 


_L3 


Fa  &aaV+Fil 


c  = 

R  \\+2m£ml-F\ 


-1 


and  C"  =  [,thk  r 


The  other  strong  realizability  condition  ( CN  >-2KR  when  F=0)  is  also  satisfied  by  this 

model.  Referred  to  as  Model-F,  the  performance  of  this  model  is  also  shown  in  Section 
5.  This  model  has  a  target  distribution  that  has  less  energy,  and  in  this  sense  it  is  similar 
to  the  simple  BGK  model  of  Section  3.  However,  unlike  the  BGK  model,  this  model  has 
the  spherical  Gaussian  as  a  solution,  is  strongly  realizable,  and  does  not  produce  a  spike 
in  the  PDF  in  the  2-component  limit.  In  addition,  unlike  the  simple  BGK  model,  the 
decay  constant,  h ,  now  depends  on  the  velocity,  v ,  and  acts  preferentially  on  the  tails  of 
the  distribution,  damping  extreme  events  more  strongly. 


5.  Model  Performance 

In  this  section  the  performance  of  these  models  is  compared  with  experimental  data  for 
return  to  isotropy.  For  each  test  case,  we  present  both  the  Reynolds  stresses  as  a  function 
of  time  and  the  Reynolds  stress  anisotropy  as  a  function  of  time.  The  anisotropy  is  the 
standard  method  for  looking  at  retum-to-isotropy,  since  it  eliminates  much  of  the 
dependence  on  the  dissipation.  However,  due  to  the  nondimensionalization  with  respect 
to  K  the  anisotropy  can  cause  errors  in  one  turbulence  component  (possibly  even 
experimental  errors)  to  appear  as  a  general  failure  of  the  entire  model.  For  this  reason  we 
retain  the  direct  Reynolds  stress  decay  plots  as  well. 

In  all  models  the  dissipation  is  determined  from  the  model  transport  equation 

f  =  -Q2-F  (C.29) 

The  value  of  Ce2  is  taken  to  be  11/6,  which  is  the  high  Reynolds  number  analytical 
solution  for  turbulence  with  a  low  wavenumber  k2  spectrum21.  In  most  of  the 
experiments  the  initial  value  of  the  dissipation  is  not  known,  and  is  obtained  by 
attempting  to  match  the  K  profile  as  well  as  possible. 

In  each  case,  we  have  solved  the  Reynolds  stress  ODE  associated  with  the  model,  using 
fourth  order  Runge-Kutta  and  very  small  time-steps.  We  have  also  solved  the 
corresponding  PDF  relaxation  models  and  obtained  very  similar  results.  However,  there 
are  further  numerical  issues  associated  with  solving  the  PDF  equations  which  we  do  not 
wish  to  address  here,  so  we  simple  present  the  ODE  results  in  this  paper. 

Because  Ce2  and  the  return  process  are  believed  to  be  Reynolds  number  dependent  we 
have  selected  only  high  Re  number  experiments  for  comparison  and  no  DNS  test  cases. 
It  must  be  noted  that  there  is  some  uncertainty  associated  with  the  experimental  results. 
First,  while  the  geometry  of  these  experiments  changes  abruptly  from  a  straining  section 
to  a  straight  section,  the  actual  cessation  of  the  mean  strain  may  not  be  quite  so  abrupt 
due  to  the  long  range  effects  of  pressure.  As  a  result,  these  decay  experiments  may  have 
some  residual  straining  in  them  at  early  times.  The  translation  of  the  zero  time  in  the  Le 
Penven  experiment,  case  III  <  0,  suggests  that  the  experimenters  were  aware  of  this 
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problem.  More  importantly,  the  initial  turbulence  for  these  experiments  has  structure, 
due  to  the  strains  imposed  to  make  the  turbulence  anisotropic.  It  is  likely  that  at  early 
times  the  relaxation  of  these  structures  also  affects  the  return  process. 

Figures  2  and  3  are  Le  Penven  et  al  cases  HI  >  0  (expansion)  and  HI  <  0  (contraction). 
Figures  4,  5  and  6  are  the  data  of  Choi  &  Lumley6  for  their  cases  A  (plane  distortion),  B 
(axisymmetric  expansion)  and  C-2  (axisymmetric  contraction)  respectively. 

Despite  the  fact  that  model-F  is  strongly  realizable  and  Model- 1  is  not,  the  two  models 
behave  very  similarly  for  all  five  experimental  test  cases.  With  the  exception  of  Figure  3 
(Le  Penven,  case  HI  <  0)  and  Figure  6  (Choi  &  Lumley,  case  C-2)  the  models  show  poor 
agreement  with  the  experimental  data,  and  tend  to  return  to  isotropy  too  quickly. 


Figure  2:  Reynolds  Stresses  and  anisotropy  for  the  case  III  >  0  from  Le  Penven ,  Gence 
and  Comte-Bellot .  Symbols  are  the  experimental  data, ,  lines  are  the  Model- 1  predictions, 
and  the  dotted  lines  are  the  Model-F  predictions. 


time  time 

Figure  3:  Reynolds  Stresses  and  anisotropy  for  case  III  <  0  from  Le  Penven,  Gence  and 
Comte-Bellot.  Symbols  are  the  experimental  data,  lines  are  the  Model- 1  predictions,  and 
the  dotted  lines  are  the  Model-F  predictions. 
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time  time 

Figure  4:  Reynolds  Stresses  and  anisotropy  for  case  A  of  Choi  and  Lumley.  Symbols  are 
the  experimental  data,  lines  are  the  Model-1  predictions,  and  the  dotted  lines  are  the 
Model-F  predictions. 
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Figure  5:  Reynolds  Stresses  and  anisotropy  for  case  B  of  Choi  and  Lumley .  Symbols  are 
the  experimental  data,  lines  are  the  Model- 1  predictions,  and  the  dotted  lines  are  the 
Mndpl-F  nrp.dicti.nnx. 
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Figure  6:  Reynolds  Stresses  and  anisotropy  for  case  C-2  of  Choi  and  Lumley.  Symbols 
are  the  experimental  data,  lines  are  the  Model- 1  predictions,  and  the  dotted  lines  are  the 
Model- F  predictions. 


6.  General  Relaxation  Models 

Rather  than  assuming  a  spherical  Gaussian  let  us  assume  that  the  anisotropic  ellipsoidal 
Gaussian  (Eq.  (C.23))  is  a  solution  to  the  relaxation  equation  (Eq.  (C.17)).  Then 

JL-— L  Wi£gB>+i5i.v'  v’  )  (C.30) 

3 1  2  J  (  det(*„m)  T  3 1  v  nv  m)  V  / 

Since  =  and  idet(/?nffl)  =  dct(Rnm)^R~J  (Jacobi’s  formula)  this 

reduces  to 

f  =  -i /(«,?' (C.31) 

Let  us  further  assume  that  -it  =  -  ( EIim  Rmj  +  E[;m  Rmi )  which  is  the  general  Reynolds 

stress  transport  model  (Eq.  (C.5)).  Then 

f =-&/(n„-nto  rj  vmv'„)  (c.32) 


This  implies  that  for  any  desired  Reynolds  stress  model,  Ely ,  a  corresponding  relaxation 
model  can  be  constructed, 


v' 


g(y)  —  CM  E[„ 


2  K 


[(2^)3det(/?nm)]1 


il/2 


h(v)  =  C, 


riffl  Rim 


M  2  K 


(C.33) 


a 

When  the  PDF  is  an  elliptic  Gaussian  we  expect  CM=  1 ,  v'„  =  v'„  =  v'„  ,  and  Rnm=  Rnm . 
The  constant  CM  can  be  a  function  of  tilde  and  hat  quantities  (such  as  up-up  and 
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up-up)  but  it  can  no  longer  be  a  function  of  the  Reynolds  stress  invariants  (like  it  was  in 

the  simpler  spherical  relaxation  model).  This  is  because  the  elliptic  Gaussian  PDF 
(unlike  the  spherical  Gaussian  PDF)  can  represent  any  state  of  the  Reynolds  stress 
invariants. 

Note  that  the  relaxation  equation  places  constraints  on  the  underlying  Reynolds  stress 
model.  It  implies  that  II(I  >  0 ,  and  II, „  R”,1  must  be  a  positive  definite  tensor. 


Conservation  of  probability  (or  mass)  is  already  satisfied  by  this  model.  Conservation  of 
momentum  requires  a  relation  between  up  and  up  (see  Appendix  L). 


(up  —«,)[:  1  +  (u„  -  un  ){um  -  um  )R: * 

=  W ^'  {  J  v ’p  v  m  V  ’«  fdy  +Rnp(Um-Um)  +  Rmp  (u„  ~  U, )} 

The  simplest  choice  is  up  =  up  then 
“p=uP+TkR^j  v'pv'mv\fdy 


(C.34) 


(C.35) 


The  choice  up  =  up  is  more  complicated  and  requires  a  symmetric  matrix  inversion 

(11,;,  RiJ  +  n,„  Rip')(un-un)  =  II,„  /£ R~l  |v >'m  v fdv .  For  certain  models  (like  the  one 

shown  below),  this  matrix  problem  is  easy  to  invert  analytically,  and  this  choice  is  also 
viable. 


The  Reynolds  stress  transport  equation  is  derived  in  Appendix  M.  Assuming  the  choice 
up  =  up  it  requires  that, 


If  the  PDF  is  an  ellipsoidal  Gaussian  then  up  =  up  (by  Eq.  (C.35)),  and  CM  =  1  (by 
definition).  In  addition,  since 

|  )  fdv  =  RmnRij  +  RnuRnj  +  RmjRni  (C.37) 

Eq.  (C.36)  gives  the  correct  limit,  Rtj  =  Ry  for  an  elliptic  Gaussian  PDF.  The  hat  and 

tilde  quantities  can  be  seen  to  be  slight  perturbations  to  the  standard  quantities  that 
precisely  account  for  any  deviation  of  the  PDF  from  an  elliptic  Gaussian  shape. 

The  model  given  by  Eqs.  (C.33)-(C.36)  represents  the  general  relaxation  model.  Using 
this  formulation,  any  Reynolds  stress  transport  model  can  also  be  implemented  as  a  PDF 
relaxation  model,  that  has  the  elliptic  Gaussian  as  a  solution.  Remember  that  =  S(j 

corresponds  to  the  case  of  no  retum-to-isotropy,  and  ny  =8ij+CR(Sij-%KRf)  is  the 
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classic  linear  retum-to-isotropy  model.  Substituting  these  expressions  into  Eqs.  (33)- 
(36)  will  produce  the  corresponding  PDF  relaxation  model.  However,  in  this  paper,  we 
do  not  wish  to  specify  II  ,  but  to  determine  what  the  general  relaxation  model  (Eqs. 

(33)-(36))  imply  about  how  it  should  be  specified. 

A 

The  general  relaxation  model  as  described  above  has  singular  h,  up,  and  Ry  in  the  two- 
component  limit  due  to  the  presence  of  Ry1 .  This  singularity  is  removed  by  the 
parameter-free  Reynolds  stress  model  II..  Making  IT,-,,  directly  proportional 

to  the  Reynolds  stress  tensor  removes  the  singularities.  The  constant  of  proportionality  is 
determined  from  the  decay  condition  TLyRy  =  2 K  (see  Section  2).  In  the  relaxation 

context  this  model  is  given  by 


e  2K2 


M  K  RmRn 


[{Iny  det(J?nm)] 


,1/2 


h  —  r  e  2g2 

K  RmRm  [2KHu-uf] 


(C.38) 


where 


(up-up)[2K+(u-u)2]  =  2Rip(ui-ui)  +  Jv'pv’(.  v',  fib  (C.39) 

Note  that  Eq.  (C.39)  is  a  particular  case  of  the  general  Eq.  (C.34)  (for  this  II(>  model).  It 

also  happens  to  be  identical  to  Eq.  (C.21),  the  general  expression  for  the  spherical 
relaxation  models  in  Section  4.  As  in  Section  4,  the  choice  of  up  =  up  or  up  =  up  is  up  to 

the  user.  For  symmetric  PDFs  if  makes  no  difference  what  the  choice  is,  since  then 
up=up=up .  For  inhomogeneous  flows,  the  PDFs  will  be  skewed  and  this  choice  may 

make  some  difference. 

For  this  model  we  also  require  the  condition  on  Ry  that, 


[4  +(«■•  -Uj)+^^][2K  +  (u-u)2] 

=  Jv  »  V  v  'jfdv  +  (un  -  un )  Jv  v )  v  'jfdv  +  Ry  (u  -  u) 


(C.40) 


This  model  (Eqs.  (C.38)-(C.40))  differs  from  those  in  Section  4  in  that  it  has  the 
ellipsoidal  Gaussian  as  a  solution. 


The  choices  for  CM  are  now  far  more  restrictive.  The  simplest  choice  is  simply  to  set 
Cu  =  1 .  Eqs.  (C.40)  and  (C.39)  are  simplified  considerably  by  choosing  up=up.  Then 


the  hat  quantities  are  defined  by  up=up+^K  jv'pV'iV)  fdy  and 


R=-L. 
lij  IK 


Jv v »  yav  u,  )(uj 
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The  equivalent  Reynolds  stress  transport  model  can  be  derived  from  this  relaxation  model 
by  assuming  the  PDF  is  an  elliptic  Gaussian.  Under  this  assumption,  the  various 
possible  choices  of  the  hat  and  tilde  quantities  are  irrelevant  and  we  find  that  all  these 
choices  are  equivalent  to, 


RisRsj 

RmnRnm 


(41) 


which  implies  the  model  parameters  are, 


c  — 4  1K 

3  R^R„, 


-land 


(42) 


We  note  that  this  model  satisfies  the  strong  realizability  constraint,  CR  - CN[-^  =  0, 

and  sits  on  the  cusp  of  the  strong  realizability  condition  CN  >  — .  In  the  2- 

component  limit,  this  model  returns  to  isotropy  as  slowly  as  physically  possible.  The 
performance  of  this  model  is  shown  below  and  it  is  referred  to  as  Model-EG  (for  elliptic 
Gaussian).  The  fact  that  the  resulting  Reynolds  stress  model  is  very  simple,  entirely 
nonlinear,  contains  no  model  parameters,  and  satisfies  strong  realizability  at  its  cusp, 
makes  Model-EG  very  intriguing. 

In  Figures  7  through  1 1  the  performance  of  Model-EG  is  compared  with  experimental 
data  for  retum-to-isotropy,  the  classic  linear  Rotta  model  (with  CR  =0.8),  and  the 
nonlinear  model  of  Sarkar  and  Speziale  ( CR  =0.7,  CN  =1.05).  The  most  interesting 
result  is  that  these  three  very  different  models  perform  very  similarly  for  all  five  test 
cases.  The  Sarkar  and  Speziale  is  slightly  better  than  the  other  two,  but  it  has  two 
adjustable  model  constants  that  were  tuned  to  exactly  these  test  cases.  The  linear  Rotta 
model  also  performs  surprisingly  well.  It  can  be  made  even  better  by  adjusting  the 
standard  value  (0.8)  downwards  (to  0.7  or  0.6).  Model-EG  matches  the  data  the  least 
well,  but  gives  quite  good  agreement  considering  there  are  no  adjustable  parameters  in 
this  model. 
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Figure  7:  Reynolds  Stresses  and  anisotropy  for  the  case  III  >  0  froth  Le  Penven,  Gence 
and  Comte-Bellot.  Symbols  are  the  experimental  data ,  lines  are  the  Rotta  Model 

nr e  dictions  /T\>  =  0  the  dotted  Imp*  nrp  the  SS  Model  nredictinns  and  dash p A  Imps  ore 


time  time 


Figure  8:  Reynolds  Stresses  and  anisotropy  for  case  III  <  0  from  Le  Penven,  Gence  and  Comte-Bellot. 
Symbols  are  the  experimental  data,  lines  are  the  Rotta  Model  predictions  (CR  =  0.8),  the  dotted  lines  are 
the  SS  Model  predictions,  and  dashed  lines  are  the  Model-EG  predictions. 
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Figure  9:  Reynolds  Stresses  and  anisotropy  for  case  A  of  Choi  and  Lumley.  Symbols  are 
the  experimental  data,  lines  are  the  Rotta  Model  predictions  fCR  =  0.8),  the  dotted  lines 
are  the  SS  Model  ore  dictions,  and  dashed  lines  are  the  ModeUEG  nredictions. 


16 


Figure  10:  Reynolds  Stresses  and  anisotropy  for  case  B  of  Choi  and  Lumley.  Symbols  are 
the  experimental  data,  lines  are  the  Rotta  Model  predictions  (CK  =  0.8),  the  dotted  lines 
nrp  thp  SS  Mndpl  nredirtinnx  nnd  dashed  linpx  nrp  thp  Mndel-F.Cr  nredirtinnx. 


Figure  11:  Reynolds  Stresses  and  anisotropy  for  case  C-2  of  Choi  and  Lumley.  Symbols 
are  the  experimental  data,  lines  are  the  Rotta  Model  predictions  (CR  =  0.8),  the  dotted 
lines  are  the  SS  Model  nredictions  nnd  dashed  lines  are  the  Model-  RG  ore  dictions. 


As  noted  earlier,  the  greatest  uncertainty  in  both  the  models  and  the  experiments  lies  in 
the  initial  conditions.  To  see  that  the  assessment  of  the  models  performance  is  not 
affected  by  these  initial  conditions  Figure  8  was  recalculated  using  a  later  time  for 
initialization.  Figure  12  shows  that  the  point  of  initialization  does  not  fundamentally 
change  the  results. 
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Figure  12:  Reynolds  Stresses  and  anisotropy  for  the  case  III  <  Ofrom  Le  Penven,  Gence 
and  Comte-Bellot,  initialized  at  0.037  seconds.  Symbols  are  the  experimental  data,  lines 
are  the  Rotta  Model  predictions  (CR  =  0.8),  the  dotted  lines  are  the  SS  Model  predictions, 
and  dashed  lines  are  the  Model-EG  predictions. 


We  conclude  this  section  by  noting  that  other  return  models  have  been  proposed  that  are 
nonlinear,  that  parameterize  CR  and  CN  as  functions  of  the  Reynolds  stress  invariants  (or 
anisotropy  invariants),  and  which  satisfy  strong  realizability6,16.  However,  these  models 
assume  that  CR  and  CN  are  polynomial  functions  of  the  invariants.  In  contrast,  the 
model  described  above  uses  linear  rational  polynomial  functions  of  the  invariants  to 
represent  the  return  parameters  CR  and  CN .  We  note  that  rational  polynomials  tend  to 
have  better  fitting  properties  than  polynomials,  and  that  the  formulated  rational 
polynomials  are  the  result  of  physical  assumptions  not  assumptions  about  functional 
behavior. 


7.  Fokker-Planck  Collision  Models 

An  alternative  to  relaxation  models  is  the  Fokker-Planck  collision  model.  This  model  is 
frequently  used  to  model  Brownian  motion,  liquid  collisions,  and  some  plasmas. 
Langevin  models  for  turbulence18'22,  are  directly  related  to  the  Fokker-Planck  equation 
and  therefore  effectively  use  this  type  of  model.  A  generalized  Fokker-Planck  collision 
operator  involves  two  as  yet  unspecified  matrices,  Gy  and  Hi} . 


9i  ~ 


a  {G9v’,f) 

+JL 

dv( 

+  av,  [ 

l  -X- 

lij  3  V, 


r 

Tt_ 


(C.43) 


The  tensor  Htj  should  be  positive  definite  for  stability  reasons.  In  Langevin  models  it  is 
convenient  to  make  Htj  isotropic  as  well.  However,  in  general  the  Fokker-Planck 

collision  model  has  considerable  flexibility  in  the  choice  of  both  the  model  tensors.  The 
model  automatically  satisfies  conservation  of  probability  and  momentum.  It  also  has  the 
ellipsoidal  Gaussian  as  a  solution23.  Multiplying  equation  (C.43)  by  v'nv'm  and 
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integrating  over  all  velocity  space  gives  the  equivalent  Reynolds  stress  transport 
equation. 

=  GmjRjn  +  GnjRjm  +  Hmn  +  H  nm  (C.44) 


By  comparing  this  with  the  generic  Reynolds  stress  transport  equation,  (Eq.  (C.5))  it  can 
be  seen  that, 


r*  .  tj  n-1  _ _ e_  T7 

Kj-  -r  n  2K  L 1  jj 


(C.45) 


In  this  way,  classic  return  models  (given  in  terms  of  Ily)  can  be  implemented  in  the 

generalized  Fokker-Planck  context.  This  transformation  is  also  discussed  in  Pope18.  The 
general  nonlinear  Reynolds  stress  return  model  (Eq.  (C.8))  is  equivalent  to 

n,  =  4  +{C„  +  C„[f -Wk]}  (4 KR:')  -  c„  4 )  (C.46) 

In  the  Fokker-Planck  context  this  implies  that 

G,  +  ^=-*[4+{c«  +  c»»-^j}(4-iJa!.j')-c»(*-lFL4)]  (C-47) 

There  are  many  possible  choices  of  Gy  and  Htj  which  satisfy  this  constraint. 


The  simplest  and  most  numerically  attractive  choice  for  HtJ  is  that  this  tensor  is 
isotropic,  Htj  =  CDeSy  where  CD  is  an  arbitrary  model  constant.  This  means  that 

Gs=-^[(X+Q+fC„)4-C„i]+f(Cs+C„[f-^]-3C0)J!8-'  (C.48) 


The  singularity  in  Gy  is  removed  by  the  particular  choice  3CD  =  CR  +  CN[j;- , 

which  is  the  choice  used  in  most  Lange vin  turbulence  models.  This  gives  the  following 
model  constants, 

ff«=[c,  +  C„(f-^)]f4  andGs=-*[(l+Ct)4+C„(f4-i)]  (C.49) 

Note  that  with  this  choice  the  realizability  constraint  CR  +CN[j--^r]  >  0  is  equivalent 
to  the  requirement  that  Hy  be  positive  definite.  Under  these  circumstances,  the  classic 
linear  return  model  (with  C„=0)  is  obtained  using  Gy  =-^(l+CR)Sy  and 
Him  =iGRSim.  Model- 1  given  in  Section  4  (with  CR  =|^-1  and  CN  =— jfe),  is 

obtained  using  Hy={-^)^dy  and  Gy  =-^(-^)[Sy  +£].  Model-F  (with 
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j"  1  RmRnm  +  p.l 
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IJ 


j^I+£na^ai_f| 


and  CV 


-i 


1 


|  ]  Rinn^nm  j? 
2  K2 


| ),  is  obtained  using  Hi}  =fSijF/[l+ -  f]  and 


Gy  =  ~2 y[^-  +  ~f]  /  [l  +  ~  F  ]  •  Note  that  in  the  2-component  limit  Htj  now  goes  to 

zero.  This  particular  splitting  (Eq.  (7.7))  will  be  unstable  in  this  limit.  Model-EG  (with 


Cr  =  3^Sr“1  and  cN=--&;>) is  obtained  usine  Hij  =°  and  GU 

This  model  is  therefore  incompatible  with  this  splitting  (unstable).  If  Hy  is  assumed  to 

be  isotropic  (and  non-zero),  then  Gy  must  become  singular  in  the  2-component  limit. 


A  more  general  splitting  is  possible  if  Htj  is  allowed  to  be  anisotropic.  Classic  Langevin 
models  require  isotropic  Hi} ,  but  the  Fokker-Planck  model  itself  only  requires  Hy  to  be 
positive  definite.  Assuming  a  positive  definite  form,  Hy  =  CDeSy  +  CE  -f  Ry  implies 


:-^(l  +  C,  +j  CN  +  2CE)SyHCs -3CD)fR:'+CN^Ry 


(C.50) 


where  Cs  =  CR  +  CN(j- ^r) ■  Again,  to  remove  the  near  singularity  CD=\CS  can  be 
chosen,  but  because  of  the  more  general  form  for  Hy  the  (realizability)  restriction 
Cs  >  0  is  no  longer  required  for  a  well  posed  model.  The  classic  linear  return  model  is 
obtained  using  Gy  =—. ykQ  +  Cr  +2CE)Sy  and  Hjm  -  f  CRSim  +CE-j?  Ry  ■  Note  that  this 
splitting  has  an  extra  free  parameter  CE  which  does  not  change  the  Reynolds  stress 
evolution,  but  does  change  the  model.  A  nonsingular  splitting  for  Model-EG  (with 
CK=i-&- 1  ^  Cn=~T&)  ™  now  given  by  Gy=-R^Ry-fCESy  and 

Hy  =fCERy ,  where  CE  is  again  an  arbitrary  parameter.  Note  that  CE  can  actually  be 
determined  by  a  dispersion  analysis  and  is  related  to  the  Kolmorgorov  constant. 


8.  Discussion 

The  retum-to-isotropy  problem  of  anisotropic  turbulence  has  been  studied  via  three  very 
different  collision  models  for  the  evolution  of  the  velocity  PDF.  The  simplest  collision 
operator  is  the  BGK  approximation  to  the  Boltzmann  collision  integral.  This  collision 
model,  -fCBGK(f-feq),  is  characterized  by  an  inverse  timescale  (which  does  not 

depend  on  the  velocity).  It  was  shown  that  if  this  model  is  to  dissipate  energy  correctly, 
the  target  state  must  have  considerably  less  energy  than  the  current  PDF  state.  Some 
models  even  use  a  target  state  with  zero  energy  (a  delta  function).  The  BGK  model 
produces  the  classic  linear  retum-to-isotropy  model,  with  the  rate  of  return  CR  =  — 

determined  by  the  energy  of  the  target  state.  The  Gaussian  PDF  is  not  a  solution  of  the 
BGK  model  even  though  theoretical  and  experimental  evidence  might  suggest  that  this  is 
desirable. 
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To  overcome  the  limitations  of  the  BGK  model,  more  general  relaxation  models  were 
constructed  in  which  the  collision  operator  g(v)-h(v)f  has  a  positive-definite  velocity 
dependent  source  term  and  a  velocity  dependent  sink  term  that  is  proportional  to  the  PDF. 
Previous  analysis  of  this  collision  model  in  the  context  of  turbulence  is  unknown  to  the 
authors.  In  Section  4  prescriptions  for  the  model  parameters  g  and  h  were  derived  such 
that  the  spherical  Gaussian  is  a  solution  to  the  evolution  equation.  Two  models  were 
derived  from  this  analysis,  Model- 1  assumed  that  the  target  distribution  has  the  same 
energy  as  the  PDF,  -|-  =  1.  It  is  only  weakly  realizable.  Model-F  assumed  that  the 

target  distribution  has  less  energy  than  the  PDF  in  the  ratio  -§-  =  F .  This  ratio  was 
chosen  because  it  makes  the  resulting  model  strongly  realizable.  While  these  models 

While  these  initial  parameter-free  relaxation  models  did  not  perform  as  well  as  might  be 
hoped,  they  set  the  stage  for  the  development  of  Model-EG.  This  model  was  shown  to 
be  the  only  nonsingular  relaxation  model  that  has  the  elliptic  Gaussian  as  a  solution.  The 
equivalent  Reynolds  stress  transport  model  is  totally  nonlinear  in  the  Reynolds  stresses 
and  was  shown  to  be  strongly  realizable.  Interestingly,  the  performance  of  Model-EG  is 
quite  similar  to  the  linear  return  to  isotropy  model.  Even  the  Sarkar  and  Speziale  model 
with  an  opposite  sign  for  the  nonlinear  term,  CN ,  performs  similarly. 

To  investigate  these  models  further,  their  trajectories  on  the  anisotropy  invariant  map 
were  plotted,  and  are  presented  in  Figure  13.  It  is  well  known  that  the  linear  Rotta  model 
has  linear  trajectories  when  plotted  on  this  anisotropy  invariant  map.  The  trajectories  of 
the  model  of  Sarkar  and  Speziale  tend  to  move  downwards  and  from  left  to  right  on  this 
map.  This  means  that  turbulence  with  two  large  Reynolds  stresses  and  one  small  stress 
will  tend  to  first  approach  a  state  with  only  one  large  stress  before  approaching  full 
isotropy.  This  implies  that  the  intermediate  stress  decays  faster  than  the  maximum  and 
minimum  stresses,  and  is  somewhat  counter  intuitive.  The  models  developed  in  this 
paper  tend  to  have  the  opposite  behavior.  Turbulence  with  one  large  stress  will  first 
decay  to  a  state  with  two  large  stresses  before  approaching  total  isotropy.  There  is  no 
experimental  data  in  the  middle  of  the  triangle  that  allows  us  to  determine  which  behavior 
is  actually  correct. 

The  top  boundary  of  the  ‘triangle’  is  the  2-component  line.  The  strongly  realizable 
models  have  trajectories  that  stay  on  this  line  and  move  to  the  left  if  they  start  on  the  2- 
component  line.  This  means  that  if  one  component  of  the  turbulence  is  zero  it  stays  zero 
for  all  time,  and  the  two  non-zero  stresses  approach  each  other  (mutual  isotropy).  This 
is  the  expected  behavior  for  two-dimensional  turbulence,  which  is  sometimes  (but  by  no 
means  always)  found  when  the  turbulence  is  two-component.  More  information  about 
the  turbulence  (than  the  Reynolds  stress)  is  clearly  necessary  to  make  return  models 
behave  correctly  in  the  2-component  limit.  Strong  realizability  seems  appropriate  when 
the  2-component  turbulence  is  also  two  dimensional,  and  weak  realizability  seems 
appropriate  otherwise. 
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Figure  13:  Invariant  triangle,  (a)  Sarkar  and  Speziale,  (b)  Model- 1,  (c)  Model-F,  (d)  Model-EG. 

Finally,  the  relationship  between  the  relaxation  models  and  the  Fokker-Planck  collision 
model  ( ~  d(C>lvi  ~  +  air  Hy  IvJ )  was  investigated.  Like  the  general  relaxation  model 

(Section  6)  the  Fokker-Planck  model  has  the  ellipsoidal  Gaussian  as  a  solution.  Because 
it  involves  derivatives  in  velocity  space,  the  Fokker-Planck  model  is  more  difficult  to 
implement  numerically  than  relaxation  models.  However,  the  Fokker-Planck  model 
(with  isotropic  Htj )  has  a  direct  correspondence  with  the  Langevin  models.  Examination 
of  Model-EG  in  this  context  showed  that  this  model  can  not  be  implemented  with 
isotropic  Htj .  Instead,  the  diffusion  coefficient  Hi}  must  be  proportional  to  Rtj . 

It  is  not  entirely  clear  that  high  Reynolds  number  decaying  anisotropic  turbulence  should 
retain  exactly  an  ellipsoidal  Gaussian  distribution  during  the  return  process.  Nevertheless 
the  assumption  of  ellipsoidal  Gaussian  evolution  is  probably  a  very  good  starting  point 
that  applies  everywhere  but  the  PDF  tails.  It  has  lead  to  the  development  of  a  novel 
parameter  free  retum-to-isotropy  model  with  interesting  properties,  and  reasonable 
agreement  with  the  data. 
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Appendix  D:  Fast  Pressure-Strain  Modeling 


The  pressure  strain  is  an  important  term  in  turbulence  modeling.  We  are  modeling  the 
two  pressure  related  terms  (pressure-strain  and  pressure-transport)  together  as  a  pressure- 
gradient  velocity  correlation  given  by  IXy  • 


_  ,  duJ)  d(P  «,;) 

dxj  dx;  dxj  dx,  _ } 

In  the  standard  approach,  FL  is  modeled  as  rapid  and  slow  pressure  strain  separately. 


<D1> 


.ar  a?-/ 

dx;  dx. 


Analysis(cf.  Chou  1945)  shows  that  in  homogenous  turbulence  the  rapid  pressure  term  in 
the  Reynolds  stress  equations  is  given  by, 


du 

IL  =2rrJLM:, 


dxq 

where  M  is  the  fourth  rank  tensor  given  by, 

M  -2  frVyjl  if  +  kpk‘  u  u)dk-  i  *»,(*)«<(*))  i  > 

M>m  ZJ  '^UiUq  +  ft*  UjUq)aK  JV  drpdrj  +  drpdr,  ’  2*lrl  aV 

The  problem  is  therefore  reduced  to  modeling  the  fourth  rank  tensor,  M,  which 
possesses  the  following  symmetry,  incompressibility,  and  amplitude  constraints, 


Mm  =M  to 
Mupq  =  0 


(D.2) 


M.  =  2  R. 

ippq  iq 

The  first  two  constraints  directly  affect  the  resulting  pressure-strain  tensor  and  must  be 
satisfied  by  any  model.  The  last  two  constraints  are  true  of  the  exact  M  tensor  but  do  not 
have  to  be  satisfied  in  order  to  obtain  a  viable  pressure-strain  model.  Durbin,  in  his  book, 
has  argued  that  satisfying  the  third  constraint  is  overly  restrictive  for  single  point  closure 
models.  Further  constraints,  such  as  realizability  of  the  pressure-strain  term  or  invariance 
to  system  rotation  in  the  2D  limit  (limited  MFI)  can  also  be  imposed.  Realizabity 
[Shuman]  requires  that  when  a  Reynolds  stress  (in  the  principal  or  diagonal  coordinate 
system)  is  zero  then  the  pressure  strain  term  does  not  cause  the  Reynolds  stress  to  be 
come  negative,  which  would  be  unphysical.  The  realizability  constraint  is  an  issue  only 
for  two  component  (2C)  turbulence  where  one  Reynolds  stress  is  zero.  It  can  be  satisfied 
by  making  certain  parameters  in  the  model  proportional  to  the  two  component  parameter 
F(which  is  0  in  the  2C  limit  and  1  in  isotropic  turbulence). 


The  standard  approach  is  to  model  M  as  a  tensor  expansion  in  terms  of  the  identity  tensor 
Sy  and  the  dimensionless  Reynolds  stress  anisotropy  tensor,  aij=^—jSiJ,  where 

K  =  \Rh  is  the  turbulent  kinetic  energy.  Launder(1975),  proposed  the  most  general 
linear  expansion.  The  anisotropy  tensor  is  not  linear  in  the  Reynolds  stresses,  but 


Ka  =  R  ~\6  K  is  linear  and  this  is  how  the  anisotropy  terms  appear  in  the  final  model. 

Higher  order  expansion  in  the  anisotropy  tensor  allows  the  model  to  satisfy  realizability 
and  all  the  incompressibility  constraints.  The  Cayley-Hamilton  theorem  allows  one  to 
truncate  the  expansion  at  third  order  in  the  tensors.  The  expansion  to  this  order  is 
performed  in  Johanson  et  al.  The  scalar  parameters  in  the  expansion  are  functions  of  the 
invariants,  and  Johannson  et  al  assume  these  functions  are  low  order  polynomials  in  the 
invariants. 

In  this  work  a  simple  linear  model  will  be  used,  to  maintain  model  simplicity.  However, 
unlike  LRR,  the  model  parameters  will  not  be  constants  and  will  be  functions  of  the 
turbulence. 


We  begin  by  looking  at  the  basic  linear  model  in  which  M  is  expanded  in  terms  of  the 
Reynolds  stress  gradients  and  identity  tensor. 

(D3) 

+  BAS„  +  B2$1R„  +  B,(RlrSj,  +  Rj,#,, ) + B,  (Rl:lSlr  +  Rustp) 

This  is  the  most  general  linear  expansion  of  M  that  satisfies  the  constraints  given  by  the 
first  constraint  in  equation  (D.2).  Imposing  the  first  incompressibility  constraint  (second 
constraint  in  equation  (D.2)  gives, 

3  A  +  2A,  +  2KB.  =  0  / 

^  ^  1  (D.4) 

3B2  +  2Z?3  +  2B4  =  0 

So  the  most  general  rapid  pressure-strain  model  in  which  M  is  a  linear  function  of  only 
the  Reynolds  stresses  contains  four  unknown  functions  and  is  given  by 

+sJ(V/, +M,  -b W+«4(V*  +v»  -l W 

If  this  model  is  going  to  be  strictly  linear  in  the  Reynolds  stress  tensor  and  the  unknown 
functions  are  only  functions  of  the  Reynolds  stress  invariants  then  the  functions  are 
highly  constrained  and  \  K  and  the  Bs  must  be  constants.  However,  in  this  model 
this  condition  will  be  relaxed  and  parameters  will  be  allowed  to  be  functions.  For 
incompressible  mean  flows  the  term  involving  5,  does  not  appear  in  the  pressure-strain 
model  and  there  are  only  three  unknown  functions. 

_  . .  ?u.  at/.  _ at/  _  at/  _ at/ 


„  du.  at/,  at/  at/ 

n-2^(a t+^)+2B^+R»X 


2&R 

3  uRPq  > 


at/, 

+2B^~^r+Rj“ 


(D.6) 


This  equation  is  can  be  presented  in  an  alternative  form, 

aj=44S#+2(BJ+B4)(V,/+V^-t'5«V»)+2(S3-S<)(V«+W)  <D-7) 

at/  at/.  du .  at/. 

where  .5,y  =^(— -£-+— -J-)and  Wtj  =^-(— — 5 - —L) .  In  general,  this  linear  model 

OXj  OXi  OXj  oxi 

satisfies  only  the  most  basic  constraints  on  the  M  tensor. 
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We  are  modeling  this  fly  as, 

n, = c;ks, + c;,(r„s„ + v*  ~iW„) + + jw 

where, 


C*isT=4A 


(D.8) 


(D.9) 


c;2=2(b3+b4) 
c;2=2(b3-b4) 

Cp2 ,  C*  and  Cwpl  are  parameters  to  be  determined. 

The  classic  model  for  the  slow  pressure  term  in  the  Reynolds  stress  equations  is  the 
Rotta(1951)  return  to  isotropy  model  given  by, 


(D.10) 


Return  to  isotropy  has  been  shown  to  occur  in  homogenous  turbulence  in  the  experiments 
of  Lumley  &  Newman  (1977).  Speziale(1991)  uses  a  C,  which  depends  on  the 
prediction.  We  follow  his  example  in  this  work  and  define, 

Ci  =  Cpl+Cp2i  (D.ll) 

£ 

where  we  found  that  Cp2  =  ~Cp2 . 

The  slow  pressure  strain  model  then  given  by, 

iC" =-<c;,  -  c;2  -\ks,)  (D.i2) 

£  A  i 

The  model  developed  above  is  for  homogenous  (or  quasi-homogenous)  turbulence.  As 
with  the  dissipation  model  it  is  likely  that  we  need  an  additional  term  near  wall  to  take 
into  consideration  inhomogeneities.  We  are  using  near  wall  term  developed  by  Amitabh 
in  his  master’s  thesis  (2001)  by  working  on  asymptotic  analysis. 

R  / 

()K  *^\1/  ) 

The  new  term  that  is  proposed  is  - v ^ - z~~ ,  adding  this  term  in  (D.12)  gives. 


dxk  dxk 


n, 


slow 


P,  £ 


dK  3(\) 


Adding  the  slow  and  fast  pressure  terms  we  get, 


(D.13) 


R  / 

,n-  n.  J* d(  Zk > 


(D.14) 


From  above  equation(D.14)  with  dissipation  tensor  (described  in  Appendix  E),  we  get 
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nr'  -*i = c',ksi + c*,2( v « + v*  -h  >  ■ +  w. + ** v 

+(l-a)^-|4)-|-fls-CV 


aVFx/aVF, 


-2i,—  Z^-2v[-^ _ s _ vs-  -  -^] 

a*>  a*>  1  f  JL "  2/r 


Combining  terms,  we  get, 

n ,"'-£,=c'/si+cyv(t«A4si)+^(W+W 


K 

+Cpl(l-a)i(|S, -f  )-f 

R  /  thULyh/E 

*  *  -M,  -  If] 
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Where, 

c;,=|(c;-c-) 

(CpI-l)(l-or)  =  C 


/>! 


(D.15) 


(D.16) 


(D.17) 

(D.18) 
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Appendix  E:  A  Model  for  the  Dissipation  Rate  Tensor  in 
Inhomogeneous  and  Anisotropic  Turbulence 

Summary 

By  including  terms  that  depend  on  gradients  a  dissipation  model  is  developed  that  is 
exact  in  the  limit  of  very  strong  inhomogeneity  (such  as  near  solid  walls  or  free-surfaces). 
Rapid  distortion  theory  (RDT)  and  equilibrium  theory  are  used  to  motivate  the 
anisotropic  terms  in  the  model.  The  resulting  model  has  only  one  free  constant  (from  the 
equilibrium  theory)  which  is  determined  via  comparisons  with  turbulent  channel  flow  at 
Re=590.  A  priori  tests  of  the  model  for  two  shear-free  boundary  layers,  channel  flow  at 
lower  Reynolds  numbers,  and  a  backward  facing  step  are  presented.  Full  simulations 
using  the  model  in  channel  flow  are  also  performed.  Comparisons  are  made  with  a 
variety  of  existing  tensor  dissipation  rate  models. 

1.  Introduction 

Reynolds  stress  transport  (RST)  equation  closures  for  turbulence  (also  referred  to  as 
single-point  second-moment  closures)  are  theoretically  capable  of  predicting  a  wide 
variety  of  complex  industrial  flows.  However,  after  many  years  of  development  RST  are 
still  not  widely  used  in  industrial  applications.  This  may  be  because  in  practice  RST 
models  often  do  not  perform  significantly  better  than  two  equation  models  in  complex 
flows.  Why  has  the  potential  of  RST  models  not  been  achieved?  One  possible 
explanation  is  that  the  development  of  RST  models  is  largely  based  on  quasi- 
homogeneous  or  quasi-isotropic  assumptions1,2.  These  assumptions  are  frequently  not 
applicable  in  engineering  flows,  particularly  those  involving  walls. 

In  this  work,  the  modeling  of  strongly  inhomogeneous  turbulence  is  explored.  In 
particular,  the  focus  of  this  paper  is  on  the  modeling  of  one  of  the  unclosed  terms  in  the 
RST  equations,  often  referred  to  as  the  dissipation  rate  tensor.  As  pointed  out  in 
Bradshaw  &  Perot ,  this  tensor  is  not  actually  equal  to  the  dissipation  rate  in 
inhomogeneous  turbulent  flows  (the  case  of  interest  in  this  paper),  so  for  brevity  and 
historical  reasons  we  simply  refer  to  this  tensor  as  the  dissipation  tensor  in  this  paper. 
Our  particular  interest  in  the  dissipation  tensor  is  due  to  the  fact  that  this  term  dominates 
in  the  region  near  a  wall.  Correct  prediction  of  the  dissipation  tensor  is  therefore  an 
important  first  step  towards  accurate  RST  models  for  complex  wall  bounded  turbulent 
flows. 


The  Reynolds  stress  transport  equation  can  be  written  as, 


DRy_ 


~(RikUj,k  +  R/tUut)+vR¥M  ^  +II„  -Tijk 


The  first  term  on  the  right  hand  side  is  the  production  term.  It  does  not  need  to  be 
modeled.  The  next  two  terms  are  the  viscous  diffusion  and  dissipation  (rate)  terms.  The 
diffusion  term  does  not  require  a  model,  and  the  dissipation  term  is  given  by 


£„  =  2vu'kUj  ~ .  This  dissipation  term  is  the  focus  of  this  paper.  The  final  two  terms,  the 
pressure- gradient  velocity  correlation  Yly  = —(p'ju' + PjUj)  >  and  the  turbulent  transport 
term  TV*  =  m'm'w'  also  require  models.  Near  a  wall,  the  turbulent  transport  is  small  and  is 

not  critical.  The  pressure-gradient  velocity  correlation  (closely  related  to  the  pressure- 
strain  term)  is  important  just  away  from  the  wall. 

Early  models  for  the  dissipation  tensor4,5  assumed  that  the  dissipation  tensor  was 
isotropic  and  given  by  the  expression  ey  =f  eSn .  Note  that  the  dissipation,  £ ,  is  a  scalar 

equal  to  one  half  of  the  trace  of  the  dissipation  tensor.  The  scalar  dissipation  is  assumed 
to  be  a  known  quantity  that  is  determined  by  its  own  transport  equation.  The  assumption 
of  isotropy  is  based  on  the  argument  that  large  velocity  derivatives  should  primarily 
occur  at  the  smallest  turbulence  scales  and  turbulence  is  thought  to  be  isotropic  at  the 
smallest  scales  (Kolmogorov6). 

While  small-scale  isotropy  of  turbulence  has  support  from  some  experiments7,  it  is 
contradicted  by  some  others810.  The  recent  theoretical  analyses  of  Hallback  et  al.n  and 
Durbin  &  Speziale12  suggest  that  under  the  action  of  mean  velocity  gradients,  even  the 
smallest  scales  and  hence  the  dissipation  tensor  must  become  anisotropic.  Brasseur13 
discusses  the  issue  in  detail. 

Since  it  is  now  widely  recognized  that  the  dissipation  tensor  is  not  isotropic  in  practice,  it 
is  often  argued  that  the  dissipation  anisotropy  should  be  modeled  along  with  the  pressure- 
gradient  velocity  correlation  following  the  practice  of  Lumley  &  Newman14.  There  is, 
indeed,  significant  evidence  to  suggest  that  the  slow  pressure-strain  correlation  and  the 
dissipation  tensor  anisotropy  are  closely  related.  However,  it  should  be  observed  that  the 
dependence  is  one  way.  The  pressure  terms  respond  to,  and  tend  to  counteract  the 
production  and  dissipation  terms.  Fast  pressure-strain  reduces  the  production  anisotropy, 
and  slow  pressure-stain  counteracts  the  dissipation  anisotropy.  In  order  to  develop 
effective  slow  pressure-strain  models  it  is  important  to  be  able  to  model  the  dissipation 
tensor  anisotropy  first. 

Some  insight  into  the  dissipation  tensor  anisotropy  in  homogeneous  turbulence  can  be 
obtained  by  using  a  Fourier  decomposition  of  the  fluctuating  velocity  field.  The 

C  <■>  ^ 

dissipation  tensor  can  then  be  written  as  etf  =  J  2 vk  u'u'j  dk.  If  the  turbulence  exists 

almost  entirely  at  one  wavenumber  magnitude  then  k2  can  be  removed  from  the  integral 
and  By  =  2vk2Ry ,  or  solving  in  terms  of  the  scalar  dissipation  ztj  =|R0 .  This  model 
was  first  proposed  by  Rotta15.  It  suggests  that  the  dissipation  anisotropy  is  equal  to  the 
Reynolds  stress  anisotropy,  ei}  =—— jSy  =x_f  =«,r  In  decaying  turbulence,  at  low 

turbulent  Reynolds  numbers  only  the  large-scale  structures  (a  single  significant  k 
magnitude)  exists  and  this  model  for  the  dissipation  tensor  becomes  exact.  The  Rotta 
model  is  therefore  frequently  referred  to  as  the  low  Reynolds  number  limit.  However,  it 
should  be  noted  that  in  many  low  turbulent  Reynolds  numbers  situations  (such  as  near 
walls)  this  critical  hypothesis  of  a  single  wave  number  magnitude  is  not  satisfied. 
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A  number  of  dissipation  tensor  models16,17  are  based  on  the  idea  of  blending  the  isotropic 
model  and  the  Rotta  model  using  a  function  that  depends  on  the  turbulent  Reynolds 
number.  These  models  have  the  form, 

£iJ  =(l-/)f  eS9  +ff  Rtj  =i£Sy  +feaiJ  (E.2) 

where  /  is  1  at  low  turbulent  Reynolds  numbers  and  0  at  high  Reynolds  numbers.  The 
model  of  Hanjalic  &  Launder17  used  /  =  l/(l  +  0.1-£) .  This  model  did  not  show  very 

good  agreement  with  DNS  data  of  channel  flow  at  Re=18018  where  the  simpler 
expression  /  =  1  (i.e.  the  Rotta  model)  was  shown  to  perform  better.  Hallback, 

19  k},2L 

Johansson  &  Burden  proposed  /  =  1/(1+^-—^)  where  Lf  is  the  integral  length  scale. 

Note  that  the  turbulent  Reynolds  number  approaches  zero  near  a  wall,  so  any  formulation 
that  uses  a  Reynolds  number  dependent  blending  function  (such  as  those  described 
above)  will  go  from  approximately  isotropic  dissipation  in  the  free-stream  to  the  Rotta 
model  near  the  wall.  An  asymptotic  expansion  of  the  dissipation  tensor  near  the  wall 
(Section  3)  shows  that  the  Rotta  model  captures  the  zeroth  order  terms  correctly  at  a  wall, 
so  these  models  show  improvement  over  pure  isotropic  dissipation  for  wall  bounded 
flows. 

Other  researchers20,21  have  proposed  using  models  other  than  the  Rotta  model  for  the  near 
wall  (or  low  Reynolds  number)  region.  These  models  have  the  form,  £tj  =f  eStj  +  fe^aU , 

where  the  wall  model  e™ll\s  trace  free.  Often,  e™11  is  defined  in  terms  of  the  wall 
normal  vector,  which  is  ill-defined  away  from  the  wall  or  at  comers.  In  addition,  in  these 
models  the  form  of  e”aU  is  formulated  specifically  for  walls  and  is  incorrect  at  a  free- 
surface  or  at  any  other  boundary  other  than  a  wall. 

While  Reynolds  number  dependent  models  capture  the  near  wall  region  better,  they  all 
revert  to  the  isotropic  model  at  high  Re  numbers  and  evidence  suggests  that  even  in  the 
high  Re  limit  the  dissipation  tensor  is  not  isotropic.  In  the  rapid  distortion  limit  Hallback 
et  al.  have  shown  that  the  dissipation  tensor  anisotropy  is  not  zero,  but  half  of  the 
Reynolds  stress  anisotropy.  The  work  of  Speziale  &  Gatski22  suggests  that  in  equilibrium 
the  dissipation  tensor  anisotropy  should  depend  on  the  shear-stress.  Finally,  Perot23  has 
shown  that  these  Reynolds  number  dependent  models  are  not  correct  for  boundaries  other 
than  walls,  such  as  slip  walls  or  free-surfaces. 

In  order  to  account  for  the  RDT  limit  Hallback,  Groth  &  Johansson11  (HGJ)  proposed  a 
nonlinear  dissipation  tensor  model.  This  model  adds  an  additional  term  proportional  to 
the  square  of  the  anisotropy  and  has  the  form 

Ey  =  f  £Sy+  fatty  +  /2*(%«„  ~  j  II Jy)  (E.3) 

where  IIa  =aijaJi  and  the  functions  are  given  by/,  =£+f  IIa  and  f2  =-f .  This  model  is 

realizable,  meaning  that  the  dissipation  tensor  in  a  certain  direction  is  zero  if  the 
turbulence  in  that  direction  is  zero.  A  similar  model  that  depends  on  the  two- 

componentality  parameter,  F  =  det(^)  was  suggested  by  Sjogren  &  Johansson24  (SJ). 
The  two-componentality  factor  F  is  1  in  isotropic  turbulence  and  0  in  two-component 
(2C)  turbulence  such  as  near  a  wall  or  a  free-surface.  In  the  SJ  model  /  =  1-0.67F  and 
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/2  =  -1.18F .  This  model  goes  to  the  Rotta  model  in  the  2C  limit,  but  does  not  satisfy 
the  RDT  condition  that  the  dissipation  anisotropy  is  half  the  Reynolds  stress  anisotropy 
under  the  action  of  large  mean  strains.  These  more  complex  models  perform  well  (away 
from  boundaries)  and  will  be  used  for  comparison  in  Section  5  where  the  model 
performance  is  evaluated. 

Speziale  &  Gatski22  have  proposed  an  algebraic  formulation  for  the  dissipation  tensor  that 
is  similar  in  construction  to  algebraic  models  for  the  Reynolds  stress  tensor.  In  the 
resulting  model  the  dissipation  tensor  anisotropy  is  solely  a  function  of  the  mean  velocity 
gradients.  Unfortunately,  the  resulting  model  reverts  to  the  (incorrect)  isotropic  model  in 
the  absence  of  mean  velocity  gradients.  This  model  is  therefore  incapable  of  representing 
the  shear-free  boundary  layers  studied  in  Section  5.  However,  the  basic  premise  of  using 
mean  velocity  gradients  to  parameterize  the  dissipation  anisotropy  (particularly  in  the 
equilibrium  limit)  is  a  reasonable  idea  which  is  adopted  later. 

Transport  equations  for  the  dissipation  tensor  can  also  be  formulated25'27.  The  Speziale  & 
Gatski  model  mentioned  above  is  a  simplification  of  such  a  transport  equation.  However, 
this  level  of  complexity  may  be  unwarranted  at  this  time  given  the  level  of  model 
uncertainty  in  the  other  RST  model  terms  (particularly  the  pressure-strain). 

In  Section  2  of  this  paper,  near  boundary  terms  for  the  dissipation  tensor  are  developed 
that  are  accurate  near  walls  and  surfaces.  These  near  wall  terms  are  derived  from  first 
principles  and  introduce  no  model  constants.  Section  3  analyzes  the  near  wall 
asymptotics  of  these  models  near  both  walls  and  free-surfaces,  and  considers  the  limit  of 
strong  inhomogeneity.  In  Section  4  the  model  development  in  regions  away  from 
boundaries  is  considered.  A  priori  tests  of  the  model  are  presented  in  Section  5  and 
compared  with  a  variety  of  existing  model  formulations.  A  brief  discussion  and 
conclusions  appear  in  Section  6. 


2.  Modeling  Strong  Inhomogeneity 

In  strongly  inhomogeneous  flows,  turbulent  correlations  such  as  the  dissipation  tensor 
change  rapidly  as  a  function  of  their  position.  Some  of  the  change  with  position  is  due  a 
change  in  the  underlying  structure  of  the  turbulence.  However,  most  of  the  change  is 
simply  due  to  the  spatial  change  in  the  turbulence  intensity.  In  the  specific  case  of  the 

dissipation  tensor,  Sy  =  2vu'kujk ,  the  dissipation  can  change  spatially  for  two  reasons. 

Either  the  gradients  correlate  differently,  or  (more  likely)  the  magnitude  of  the  velocity 
fluctuations  has  simply  changed.  These  different  effects  can  be  isolated  by  using  the 
following  decomposition.  Let  the  fluctuating  velocity  be  decomposed  as  u'  =  Q.Ju. .  The 

tensor  Qyis  assumed  to  be  a  known  quantity  (related  to  the  velocity  fluctuation 

magnitude).  It  is  an  average  quantity  and  does  not  change  in  time  for  statistically  steady 
flows  or  along  homogeneous  directions.  The  underlying  temporal  and  spatial  fluctuations 
of  the  velocity  field  are  captured  by  the  dimensionless  quantity  u,.  Changes  in  the 
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dissipation  due  to  changes  in  the  turbulence  magnitude  will  be  captured  by  Qy .  Changes 
in  the  underlying  turbulent  structure  will  be  manifest  in  u. . 


Substituting  this  formula  into  the  equation  for  the  Reynolds  stress  tensor  gives  a 
relationship  between  the  structure  correlation  and  the  Reynolds  stress  tensor. 


Rn  =  u'u':  =  OJL0JL  =  Q,„umJO , 

ij  i  j  zz'in  nxZ'jm  m  *-,in  n  m^jm 


(E.4) 


The  magnitude  tensor  is  not  a  fluctuating  quantity  and  therefore  can  come  out  of  the 
average.  Substituting  this  decomposition  info  the  dissipation  tensor  formula  gives, 


eij  =  2VUi,k“j,k  =  2HQinA+Q<n“n,k)(Q>m,k“n,+QjmUm,k) 

~  2v(QinJ,Q  jmk  \  (QinQ  jm  )*  (unfin),k 


~^QinQjm^n,k^m,k  ~^"2^Qin,kQjm  QinQ  jm,kQ^m,kMn  ,jfc )  } 


(E.5) 


If  it  is  required  that  Qy  be  invertible  then  the  first  two  terms  in  the  expression  can  be 


found  from  Eqn.  (E.4)  and  are  exact.  The  third  term  is  the  dissipation  of  the  velocity 
structure.  It  requires  a  model.  However,  the  velocity  structure  is  quasi-homogeneous 
(by  design),  and  so  standard  dissipation  models  are  expected  to  perform  well  in  this 
context.  The  final  term  is  the  product  of  two  differences.  It  is  assumed  to  be  small  and 
evidence  to  that  effect  can  be  found  in  Perot  &  Moin.28  In  regions  of  strong 
inhomogeneity  the  first  term  dominates  and  Eqn.  (E.5)  becomes  exact  irrespective  of  the 
model  used  for  the  third  term  in  Eqn.  (E.5)  or  the  size  of  the  fourth  term. 


One  possible  definition  for  Q..  is  that  it  represents  all  the  magnitude  information  (Perot, 

&  Moin,29).  In  this  case  unum  =  5nm  and  Eqn  (E.4)  becomes  Ry  =  QinQjn  or  R  =  QQT  and 

Q  is  the  matrix  square  root  of  the  Reynolds  stress  tensor.  This  definition  of  Q  is  actually 
not  unique,  Q  can  be  symmetric  or  lower  triangular  for  example.  The  symmetric  square 
root  however,  seems  to  be  the  most  natural.  Like  regular  square  roots,  the  sign  of  Q  is 
also  not  well  defined.  Since  Q  always  appears  in  pairs,  this  distinction  is  not  important. 
With  this  definition  of  Q  ,  the  second  term  of  Eqn  (E.5)  is  identically  zero,  and  the  model 
is  given  by 

Eij  =  2VQin,kQjn,k  +  QinKmQjm  (E.6) 

where  the  fourth  term  of  Eqn.  (E.5)  is  assumed  to  be  negligible. 


This  near  wall  model  is  elegant,  but  inconvenient  to  implement.  Finding  Q  requires 
determining  the  eigenvectors  and  eigenvalues  of  R.  In  this  paper  we  consider  a  simpler 
implementation  of  Eqn  (E.5).  In  order  to  gain  implementation  simplicity  we  therefore 
assume  that  Q  is  isotropic  and  is  scaled  by  the  turbulent  kinetic  energy,  Q..  =  KmSy . 

With  this  definition  of  Q,  Eqn  (E.4)  gives  the  relation  =  and  the  fourth  term  in 
Eqn.  (E.5)  is  identically  zero.  Eqn.  (E.5)  then  becomes, 

t„=2y(K'l\(,K''\^+vKJ^),+Ki,  (E.7a) 

which  is  an  exact  relation.  This  can  alternatively  be  written  as 

(E.7b) 


e,y=2 v{K 


1/2 ) 

)’n  \  )  „ 


+  K~Eij 
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Note  that  Eqn.  (E.7)  only  becomes  a  dissipation  tensor  model  when  a  quasi-homogeneous 
dissipation  tensor  model  (for  Ke^ )  is  hypothesized.  The  quasi-homogeneous  dissipation 

tensor  should  be  significantly  easier  to  model  than  the  dissipation  tensor  itself.  The 
quasi-homogeneous  model  is  discussed  in  section  4.  In  the  next  section,  the  near  wall 
behavior  of  Eqn.  (E.7)  is  analyzed. 


3.  Asymptotic  Analysis  Near  Boundaries 

The  behavior  of  turbulence  quantities  near  a  boundary  can  be  determined  by  using  Taylor 
series  expansions  in  the  coordinate  direction  normal  to  the  boundary  (Launder  & 
Reynolds30).  Using  the  convention  that  y  is  the  direction  normal  to  a  wall  the  fluctuating 
velocity  can  be  expanded  as, 

u'  =  ai{x,z,t)  +  ybi(x,z,t)  +  y2Ci(x,z,t)  +  ...  (E.8) 

At  a  solid  wall  the  velocity  goes  to  zero,  so  all  the  at  are  zero.  Continuity  applied  very 
close  to  the  wall  implies  b2=  0 . 


Substituting  Eqn  (E.8)  into  the  definition  for  the  dissipation  tells  us  that  near  a  wall, 
e,  ,  =  2v{  (ft,2 )  +  y(4£,c,  )+y2  (4c,2  +  6  bldl  +bul2  +£>132)+...} 

el2  =  2v{y{7bxci)+  y2  (3fe,42  +  4c,c2)  +  y3  (4bte2  +  6c,42  +  64,c2  +fcuc21  +fc13c23)+...} 


e22  =  2v{  y2  (4c22  )  +  y3  (1 2 c2d2 )  +  y4  (9 d2  + 1 6 c2e2  +  c2  2  +  c2  32)  + . . . } 

The  £33  component  behaves  just  like  £n.  A  similar  expansion  for  the  Reynolds  stress 
tensor  can  also  be  performed. 

Rll  =  y\b2)  +  y3(2blcl)+... 

Rn  =  y\blc2)  +  y\bld1+clc1)+...  (E.  10) 

^22  =  y*(c22)  +  ys(2c2d2)+... 

The  leading  order  terms  in  the  dissipation  tensor  and  the  Reynolds  stress  tensor  are  very 
similar.  However,  the  coefficient  is  different  in  each  case.  Rotta’s  model  gets  the  0(1) 
terms  of  the  dissipation  tensor  correctly  (i.e.  the  leading  term  of  the  eu  and  e33 

expansion),  but  it  will  under  predict  the  leading  order  terms  of  the  other  two  dissipation 
components.  Although  the  wall  is  at  a  low  turbulent  Reynolds  number,  Rotta’s  model 
does  not  work  entirely  correctly.  The  amplitude  of  fluctuations  normal  to  the  wall  and 
those  parallel  to  the  wall  are  very  different,  and  the  basic  assumptions  used  to  derive  the 
Rotta  model  are  violated. 


Even  if  leading  order  terms  of  en  and  £22  are  wrong,  does  it  matter?  They  still  go  to 
zero  at  the  wall.  Interestingly,  if  wall  functions  are  not  used  it  matters  a  great  deal  (using 
wall  functions  with  a  RST  model  largely  defeats  the  purpose  of  having  a  RST  model,  see 
Speziale31).  Near  the  wall,  the  dissipation  and  pressure-gradient  velocity  correlation 
exactly  balance  the  diffusion  term.  If  the  leading  order  behavior  of  the  dissipation  is 
incorrect,  the  Reynolds  stresses  are  too  large  near  the  wall  and  as  a  result  they  are  also 
too  large  away  from  the  wall.  Trying  to  reduce  these  Reynolds  stresses  via  terms  in  the 
model  (rather  than  fixing  the  root  cause)  often  leads  to  instability  in  the  wall  bounded 
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RST  equation  system.  Note  that  one  reason  elliptic  relaxation  models  work  well  has 
nothing  to  do  with  ellipticity.  These  models  allow  an  extra  boundary  condition  to  be 
imposed  (because  they  hypothesize  an  extra  equation).  This  additional  boundary 
condition  forces  the  correct  near  wall  behavior  of  the  Reynolds  stresses.  In  essence,  the 
elliptic  relation  forces  the  near  wall  behavior  of  the  dissipation  tensor  to  be  correct  via 
additional  boundary  conditions.  In  standard  RST  models  (where  six  additional  equations 
and  their  boundary  conditions  are  not  available),  correct  leading  order  behavior  of  each 
dissipation  term  is  highly  desirable. 


As  mentioned  earlier,  it  is  also  possible  to  formulate  models  with  the  correct  near  wall 
asymptotics  by  using  the  wall  normal  vector  or  distance  to  the  wall  along  with  a  blending 
function.  This  works,  and  is  standard  practice,  but  these  models  have  serious  deficiencies 
in  their  generality.  Typically  they  work  only  at  walls. 


The  boundary  conditions  at  a  slip  wall  (or  stationary  free-surface),  impose  different 
constraints  on  the  expansion.  We  now  find  that  bx  =  a2  =  b3  =  0 ,  and  continuity  implies 

Ch.x+a3,z+b2=:^‘ 

At  a  stationary  free-surface  the  dissipation  behaves  as, 
e„  =2v{(au2  +al2)  +  y2(axxcxx  +a13c,3  +4c,2)+...} 

Ei2  =  2v{  y(aub2A  +  <3,3^3  +  2  cxb2 ) + y2  (a,  ,c2J  +  al3c2i  +  3djb2  +  4qc~2) + ... }  (E.l  1) 

E21  =  2v{(b2)  +  y(4c2b2)+...} 
and  the  Reynolds  stress  tensor  is, 

Rll=(a2)  +  y2(2aicl)+... 

Hi2  =  y(.ail>2)+y2(alc2)+...  (E.12) 

R22  =  y2(b22)+y\2b2c2)+... 

At  a  free-surface  there  is  no  longer  a  clear  relationship  between  the  dissipation  tensor  and 
the  Reynolds  stress  tensor.  Rotta’s  model  will  cause  e22  to  be  zero  at  the  surface  when  it 
should  be  finite.  Also  note  that  a  free-surface  is  no  longer  a  low  turbulent  Reynolds 
number  situation,  so  blending  models  (Eqn  E.2)  will  produce  the  isotropic  limit  near  the 
surface.  The  isotropic  model  does  give  a  finite  value  for  e22  but  it  will  be  shown  in 

Section  5  that  it  is  far  too  large,  and  that  the  dissipation  near  a  free-surface  is  not  close  to 
isotropic. 


The  near  boundary  behavior  of  the  proposed  model  can  be  determined  from  the  behavior 
of  the  Reynolds  stresses.  For  a  no-slip  wall  we  find  that 

K  =  y2  2  (V  + ■ h2)  +  y3  j  Wa  +  2 b&) + ...  (E.  1 3) 

and 


4{l+yl^=j!C3)+. ..} 

y  (V+ij2) 


plugging  into  the  model  equation  (Eqn.  E.7)  gives, 


(E.14) 
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(E.15) 


£n  =2v((bl2)  +  y(4blcl)  +  0(y*)}  +  K£n 
el2  =  2v{  y(2 b^) + y 2  (3 +3  7^+b^  +  0(y3) } +  Ken 

e22  =  2v(y^)  +  y\^2+2^^)  +  0(/)}  +  Ki22 

So  the  proposed  expression  for  the  dissipation  tensor  (Eqn  E.7)  captures  the  0(1)  and 
0(y)  terms  exactly  and  at  least  75%  of  the  0(y2)  terms,  when  implemented  near  a  wall. 
Since  K  is  0(y2)  this  analysis  shows  that  the  quasi-homogeneous  dissipation  model  can 
be  as  high  as  0(1)  at  the  walls,  without  affecting  the  near  wall  asymptotics  described 
above.  Before  considering  the  behavior  of  the  quasi-homogeneous  dissipation  tensor  in 
any  more  detail,  let  us  consider  the  behavior  of  the  proposed  decomposition  (Eqn.  E.7) 
near  a  free  surface. 

Near  a  free  surface  the  kinetic  energy  is  given  by 

K  =  +  a*) +  y2 \{b2c2  +  2axcx  +2  a3c3)+...  (E.16) 

and  _ 

(If)2  ={[i(?+?)(1]2+[i(?+?)>3]2}/[(^2  +  «32)]2  +  ^(/)  (E.17a) 

we  can  also  show  that 

K,&h=0(l) 

Kk(^\t=0(y)  (E.17b) 

Kk(^\k=0(y2) 

so  the  near  boundary  terms  in  Eqn.  (E.7)  have  the  same  type  of  behavior.  This  requires 
the  £12  model  to  go  like  O(y)  near  the  surface  and  £22  to  be  0(1).  Looking  at  the  exact 
expressions  for  the  dissipation  tensor  near  a  free-surface  it  is  clear  that  capturing  the 
leading  order  £n  and  £n  terms  exactly  is  not  possible.  Derivative  information  is  not 
available  to  a  RST  model.  However,  the  leading  two  terms  of  the  £22  expression  can,  in 
theory,  be  obtained  exactly  from  Reynolds  stress  information.  Also  note  that  the  error  in 
both  the  wall  and  free-surface  expressions  for  £22  can  be  represented  by  ■^■R22 .  The  error 

is  the  same  irrespective  of  the  boundary  type.  At  both  boundaries,  the  flow  becomes  two- 
component,  so  we  will  use  the  2C  parameter  F  to  model  this  missing  contribution  for  £22 . 
This  extra  term  is  2v(F'n) n(F112) mRnmSiJ / F .  Technically  we  are  now  modeling  the  quasi- 
homogeneous  dissipation,  Key .  This  is  the  near  boundary  contribution  of  the  quasi- 

homogeneous  dissipation  due  to  the  2C  nature  of  the  turbulence  near  these  boundaries. 
This  term  is  higher  order  for  £u  for  £n  terms  near  both  walls  and  free-surfaces,  and  so  it 
only  affects  the  £22  dissipation  component.  At  a  solid  wall,  this  enhancement  has  only  a 
very  small  affect  on  the  model.  However,  at  a  free  surface  the  2C  affects  can  be  seen 
very  clearly.  The  importance  of  this  2C  correction  is  demonstrated  in  Section  5. 
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4.  Quasi-Homogeneous  Dissipation 

In  homogeneous  turbulence,  the  boundary  (or  gradient)  terms  drop  out  entirely  and  the 
quasi-homogeneous  dissipation  remains  to  be  modeled.  Hallback  et  al.19  show  that  in 
initially  isotropic  homogeneous  turbulence  the  dissipation  anisotropy  should  be  half  of 
the  Reynolds  stress  anisotropy  under  the  action  of  rapid  irrotational  strain  or  shear.  This 
will  be  referred  to  as  the  RDT  limit.  The  experiments  of  Crow32  and  Lee  &  Reynolds33 
show  that  this  ratio  does  not  remain  Vi  when  the  turbulence  is  anisotropic,  and  in  the 
extreme  limit  of  axisymmetric  2C  turbulence  it  is  seen  to  be  close  to  1  (which  is  the  Rotta 
model). 

The  practice  of  expanding  model  parameters  in  polynomial  expansions  of  the  potential 
unknowns  is  a  rational  way  to  proceed,  and  is  certainly  viable  when  the  unknowns  are 
known  to  be  small.  However,  when  the  objective  is  to  capture  an  entire  functional  range 
the  use  of  polynomial  expansions  can  be  detrimental.  Rational  polynomials  have  a 
greater  fitting  capability.  In  this  work,  we  propose  a  simple  tensorally  linear  model  for 
the  quasi-homogeneous  dissipation,  in  which  the  blending  parameter/ is  a  function  of  F. 
This  is  similar  to  the  models  of  Johansson,  however  we  hypothesize  a  rational  polynomial 
expansion,  f  =-&?,  rather  than  a  simple  polynomial  series.  This  results  in  the  quasi- 

homogeneous  model  eij=e(^Sij+1^^-)  =  £^Sij+1^aiJ)  and  In 

isotropic  turbulence,  this  model  gives  the  correct  RDT  anisotropy  ratio  of  Vi,  and  in  2C 
turbulence  it  gives  the  correct  anisotropy  ratio  of  1.  In  theory,  a  slightly  more  complex 
blending  might  be  desired  in  which  /  =  )F  where  the  function  g  goes  to  zero  as  the 

turbulent  Reynolds  number  becomes  small  and  approaches  1  at  high  turbulent  Reynolds 
numbers.  We  have  not  pursued  this  added  level  of  complexity  at  this  time. 

Finally,  we  note  that  the  dissipation  anisotropy  could  be  a  function  of  the  mean  flow 
gradients,  not  just  the  Reynolds  stress  anisotropy.  Typically,  dissipation  anisotropy  is 
not  modeled  in  this  way  because  one  does  not  expect  sudden  changes  in  the  mean  flow  to 
have  an  instantaneous  affect  of  the  dissipation.  However,  in  equilibrium  situations,  there 
could  be  a  good  correlation  between  the  two  tensors.  Reynolds  stresses  are  frequently 
modeled  using  this  type  of  hypothesis  (eddy  viscosity  hypothesis  of  Bousinesq).  In  fact, 
models  which  only  depend  on  the  Reynolds  stress  anisotropy  will  have  the  dissipation 
anisotropy  aligned  along  the  same  principal  directions  as  the  Reynolds  stress  anisotropy. 
We  know  that  these  anisotropy  directions  are  not  always  aligned  (in  channel  flow  they 
disagree  by  8  degrees  at  y+  =30).  In  this  work  we  therefore  hypothesize  that  the  quasi- 
homogeneous  dissipation  tensor  can  also  be  a  linear  function  of  the  shear-stress  tensor, 


The  model  for  the  dissipation  tensor  then  becomes, 


Sy+C'KS,!  (E.18) 


^=2v(r-,)2R(y+v^(f)it+£^f^X+f~^^  +  2v- 
The  single  parameterC*  =0.18F/(1+F)2  is  set  by  comparing  the  £12  component  of 
turbulent  channel  flow  at  Re=590.  In  theory,  the  constant  C*  should  be  a  function  of 
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jif ,  such  that  C*  goes  to  zero  when  -§g  is  zero  (the  RDT  limit).  We  have  not  explored 
this  level  of  detail  in  this  work. 

The  scalar  e  =  \en  is  the  trace  of  the  quasi-homogeneous  dissipation.  It  has  units  of 
inverse  time  or  frequency  and  can  be  obtained  by  taking  one  half  the  trace  of  Eqn.  1 8, 
Ke-e- 2 v(KU2,„ )2 - 3v .  The  quasi-homogeneous  dissipation  £  (or  its 

closely  related  form  i  =  Ke)  is  an  interesting  inverse  timescale  that  has  been  used 
previously  in  some  near  wall  turbulence  models  (e.g.  Launder  &  Shima,32).  It  is 
attractive  because  at  a  wall  it  is  finite,  whereas  the  standard  inverse  timescale  f  is 
singular  and  goes  like  y2  at  a  wall.  Note  that  from  its  definition,  e  is  a  positive  quantity. 
However,  due  to  numerical  inaccuracy  in  the  calculation  of  gradients,  calculating  e  from 
the  formula  above  can  lead  to  large  errors  or  negative  values  when  implemented  on  a 
*  computer.  In  practical  implementation  either  a  transport  equation  is  solved  directly  for  £ 
rather  than  the  more  common  £  transport  equation  (as  in  many  low  Re  number  k/e 
models),  or  we  sometimes  use  g  =  f(1+IO(<|V(1A.i/2)l/A:)  to  guarantee  a  positive  inverse 

timescale  with  finite  near  wall  behavior. 

While  the  proposed  model  (Eqn.  E.18)  looks  somewhat  complex,  it  is  relatively  easy  to 
implement.  Many  of  the  terms  combine  with  similar  looking  terms  in  the  pressure-strain 
model,  and  if  the  Reynolds  stress  anisotropy  equation  is  solved  rather  than  the  RST 
equation,  then  some  of  the  terms  drop  out  or  simplify  even  further. 

It  is  important  when  implementing  this  model  to  have  a  numerical  method  that  is  capable 
of  accurately  calculating  gradients.  At  the  wall,  certain  terms  should  exactly  balance. 
Numerically  they  will  only  approximately  balance  and  if  the  disagreement  is  large 
enough,  the  numerical  implementation  (not  the  model)  becomes  unstable.  Quantities 
with  high  power  law  behavior  (R22  -0(y4))  can  be  quite  hard  to  differentiate  accurately 
with  low  order  numerical  methods.  For  this  reason,  the  anisotropy  equations  (rather  than 
Reynolds  stress  equations)  are  somewhat  easier  to  solve  with  low  order  numerical 
methods. 


5.  Model  Results 

In  this  section,  the  proposed  dissipation  model  (Eqn.  E.18)  is  compared  against 
experimental  and  DNS  data.  The  performance  of  the  model  is  compared  to  a  number  of 
other  dissipation  tensor  models  that  have  been  mentioned  in  the  text.  The  majority  of  the 
tests  are  a  priori  tests  using  data  for  the  Reynolds  stresses  and  dissipation  plugged 
directly  into  Eqn  18. 

These  tests  are  a  useful  way  to  directly  isolate  if  the  model  can  represent  the  dissipation 
tensor  accurately.  However,  it  is  possible  to  construct  models,  which  perform  well  in  a 
priori  tests  but  do  not  perform  well  in  practice.  These  models  are  unstable  and  move 
away  from  the  desired  solution  rather  than  towards  it.  To  demonstrate  stability  we  will 
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also  present  at  the  end  of  this  section  some  solutions  of  turbulent  channel  flow  that  use 
this  dissipation  tensor  model  in  a  full  RST  prediction. 

Our  first  test  case  does  not  involve  the  inhomogeneous  terms  at  all.  It  is  a  test  of  the 
quasi-homogeneous  part  of  the  model.  Figure  1  shows  the  model  performance  in 
axisymmetric  rapid  contraction  of  homogeneous  turbulence.  In  this  flow  the  turbulence 
is  initially  isotropic  and  becomes  increasingly  2C  as  time  proceeds.  Because  the 
turbulence  is  axisymmetric  only  one  component  of  the  dissipation  need  be  analyzed.  The 
figure  shows  the  en  component  as  a  function  of  a,, ,  at  various  times  during  the 
simulation  (experiment).  The  isotropic  model  gives  a  flat  line,  and  the  dashed  line  with  a 
slope  of  1  is  the  Rotta  model.  This  is  a  relatively  high  Reynolds  number  test  case,  so 
models  that  switch  between  the  isotropic  model  and  the  Rotta  model  based  on  a  blending 
parameter  which  is  a  function  of  the  turbulent  Reynolds  number  (most  models)  will  be 
essentially  isotropic  (very  close  to  a  horizontal  line  through  the  origin).  The  HGJ  model 
was  designed  for  this  flow  and  therefore  performs  well  for  this  case.  The  SJ  model  (not 
shown)  performs  very  like  HGJ. 


Figure  1.  Dissipation  anisotropy  in  axisymmetric  contraction.  Open  triangles 
denote  the  experimental  data  of  Crow32  ( S*2(t  =  0)  =  0.5  —  2.0,  Re^  =  15  —  100) , 
open  circles  denote  the  experimental  data  of  Lee  &  Reynolds33 
( $2  =  0)  **  0-97  to.0.71,ReA  «50) ,  thick  dashed  line  denotes  HGJ  model,  thin 

line  denotes  isotropic  model,  thin  dashed  line  denotes  Rotta  model,  and  thick  line 
denotes  the  proposed  model. 

Next,  we  wish  to  examine  the  inhomogeneous  terms.  The  quasi-homogeneous  term 
cannot  be  completely  eliminated  in  any  flow,  but  shear-free  boundary  layers  provide  an 
opportunity  to  evaluate  the  model  in  a  strongly  inhomogeneous  situation  with  few  other 
complicating  effects.  In  Figure  2,  the  model  predictions  of  a  shear-free  boundary  layer 
next  to  a  solid  wall  are  compared  with  DNS  data35.  Like  the  previous  case,  this  flow  is 
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axisymmetric  and  time  developing,  but  unlike  the  previous  case,  it  also  has  strong 
gradients  in  the  direction  normal  to  the  wall.  The  figure  shows  the  e22  dissipation  at  two 
different  times  after  the  wall  appears.  As  predicted  by  the  asymptotic  analysis,  the  Rotta 
model  is  too  small  near  the  wall.  The  HGJ  model  transitions  from  a  combination  of  Vi 
isotropic  and  Vi  Rotta  well  away  from  the  wall  to  all  Rotta  near  the  wall  where  the 
turbulence  is  2C.  The  SJ  model  (not  shown)  is  almost  identical  to  HGJ.  Other  models 
which  blend  the  isotropic  and  Rotta  models  based  on  the  Reynolds  number  will  behave 
like  the  Rotta  model  near  the  wall.  All  but  the  proposed  model,  underpredict  the  normal 
dissipation  component  near  the  wall. 


y 

a)  After  time  =  2.01 
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b)  After  time  =  6.04 


Figure  2.  Shear-free  turbulent  boundary  layer  next  to  solid  wall.  Circles  denote 
DNS  data  of  Perot  &  Moin35,  thick  dashed  line  denotes  isotropic  model,  thin  dashed 
line  denotes  Rotta  model,  thin  line  denotes  HGJ  model  and  thick  line  denotes 
proposed  model. 


The  shear-free  boundary  layer  next  to  stationary  free  surface  is  shown  in  Figure  3.  Again 
two  times  are  shown,  and  the  definitions  of  the  lines  and  symbols  are  the  same  as  in 
Figure  2.  This  flow  is  no  longer  low  Reynolds  number  near  the  surface,  and  so  the 
under-prediction  of  the  Rotta  and  HGJ  models  is  even  more  obvious  in  this  case.  Most 
other  models  will  behave  like  the  isotropic  model  for  this  flow.  The  proposed  model 
captures  the  near-surface  dissipation  correctly  using  no  adjustable  constants. 

In  figure  4,  the  various  terms  in  the  present  model  are  split  out  so  the  magnitude  and 
location  of  each  contribution  can  be  ascertained.  Figure  4(a)  is  the  shear-free  surface  at 
time  2.0  and  figure  4(b)  is  the  shear-free  wall  at  time  2.0.  It  is  clear  that  the  term 

involving  2v — ^-jr — -  is  critical  in  the  free-surface  case. 

In  Figure  5  the  model  is  tested  in  turbulent  channel  flow  at  Re  59036.  The  proposed 
model  performs  well  for  all  the  dissipation  components.  The  other  models  have  difficulty 
predicting  thef12  and  e22  components.  The  value  of  C *  was  tuned  for  en  in  this  case. 
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Figure  3.  Shear-free  turbulent  boundary  layer  next  to  free-surface.  (See  Figure  2 
for  caption). 
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a)  Shear-free  surface 
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b)  Shear-free  wall 


Figure  4.  Each  term  in  equation  18  for  shear-free  turbulent  boundary  layers  at 
time  2.0.  The  chain  dotted  line  denotes  2v{Kvl,nfRij  +  £-^FjSijK+i-^FRiJ ,  thin 

line  denotes  vKk  (^)  ,  and  thick  line  denotes  2v^^f^LKSij .  The  shear-stress 

'  '  ,Jk  '  * 

term  is  zero  for  both  these  flows. 
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Figure  5.  Dissipation  tensor  in  channel  flow  (Re  =  590).  Circles  denote  DNS  data 
from  Ref  36,  long  dashed  line  denotes  SJ  model,  chained  dashed  line  denotes  HGJ 
model  and  thick  line  denotes  proposed  model.  (In  the  figure  of  en ,  thick  dashed  line 
represents  the  proposed  model  without  the  shear  stress  term) 


With  C*  =  0 ,  the  small  dashed  line  is  obtained.  The  model  predictions  for  channel  flow 
at  Re=395  and  Re=180  are  shown  in  Figure  6  and  Figure  7  respectively.  The  lowest 
Reynolds  number  case  shows  some  discrepancies  for  £n .  This  might  be  fixed  by  making 

C*  a  function  of  e/SK  as  suggested  earlier. 

The  behavior  of  the  model  in  rotating  channel  flow  is  shown  in  figure  8.  The  test  case  is 
the  Ro=0.15  DNS  case  of  Andersson  and  R.  Kristoffersen37  at  a  Re  of  194.  The  model 
does  a  good  job  of  predicting  the  very  different  behaviors  on  each  side  of  the  channel 
(especially  given  the  very  low  Re  number  of  the  simulation). 

The  proposed  model  is  compared  with  DNS  data  for  the  flow  past  a  backward  facing  step 
in  Figures  9,  10  and  11.  Several  locations  are  shown,  both  before  (4h)  and  near  (6h) 
reattachment,  and  well  downstream  (lOh)  during  the  boundary  layer  recovery.  The 
figures  show  good  agreement  with  the  DNS  data38  in  the  turning  mixing  layer.  The 
boundary  layer  near  the  wall  is  more  difficult  to  see,  but  behaves  similarly  to  the  previous 
channel  flow  results.  In  the  mixing  layer,  this  tensoraly  linear  model  performs  similarly 
(or  better  in  the  en  case)  to  the  more  complex  nonlinear  model  of  HGJ  and  is  much  more 
accurate  than  the  isotropic  and  Rotta  models. 
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Figure  6.  Dissipation  anisotropies  in  channel  Flow  (Re  =  395).  Circles  denote 
DNS  data  from  Ref  36  and  thick  line  denotes  nronosed  model. 


It  is  interesting  to  quantify  the  effect  of  the  nonlinear  term  in  the  SJ  and  HGJ  models.  In 
figure  12,  these  models  are  shown  with  f2  set  to  zero,  for  the  7h  downstream  location  on 
the  backward-facing  step.  The  contribution  of  the  nonlinear  term  is  not  that  large,  given 


Figure  7.  Dissipation  anisotropies  in  channel  Flow  (Re  =  180).  Circles  denote 
DNS  data  from  Ref  36  and  thick  line  denotes  proposed  model. 
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Figure  8.  Dissipation  anisotropies  in  rotating  channel  Flow  (Ro=0.15,  Re  =  194). 
Circles  denote  DNS  data  from  Ref  37  and  thick  line  denotes  proposed  model.  Thin 
line  in  the  last  figure  is  the  model  with  C*=0. 
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Figure  9.  Dissipation  tensor  in  the  flow  over  backward-facing  step  (at  4h). 
Circles  denote  DNS  data  from  Ref  38,  long  dashed  line  denotes  the  isotropic 
model,  chain-dotted  line  denotes  the  Rotta  model,  thin  line  denotes  the  SJ  model, 
the  small  dashed  line  denotes  the  HGJ  model,  and  the  thick  line  denotes  the 
proposed  model. 


Figure  10.  Dissipation  tensor  in  the  flow  over  backward-facing  step  (at  6h). 
Circles  denote  DNS  data  from  Ref  38,  long  dashed  line  denotes  the  isotropic 
model,  chain-dotted  line  denotes  the  Rotta  model,  thin  line  denotes  the  SJ  model, 
the  small  dashed  line  denotes  the  HGJ  model,  and  the  thick  line  denotes  the 
proposed  model. 


Figure  11.  Dissipation  tensor  in  the  flow  over  backward-facing  step  (at  lOh). 
Circles  denote  DNS  data  from  Ref  38,  long  dashed  line  denotes  the  isotropic 
model,  chain-dotted  line  denotes  the  Rotta  model,  thin  line  denotes  the  SJ  model, 
the  small  dashed  line  denotes  the  HGJ  model,  and  the  thick  line  denotes  the 
proposed  model. 
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Figure  12.  Dissipation  anisotropies  in  flow  over  backward-facing  step  (at  6h). 
Circles  denote  DNS  data  from  Ref.  38,  chained  dashed  line  denotes  SJ  model, 
dashed  line  denotes  HGJ  model,  thin  line  denotes  S J  model  without  nonlinear  term, 
thick  line  denotes  HGJ  model  without  nonlinear  term. 


its  added  complexity  we  have  chosen  not  to  include  such  a  nonlinear  term  in  the  proposed 
model.  The  choice  of  the  function  f  =  allows  us  to  satisfy  the  RDT  limit  and 
realizability  without  the  nonlinear  term  present. 

A  priori  tests  such  as  those  described  above  can  be  very  informative  about  the  quality  of  a 
model.  Nevertheless,  it  is  possible  to  formulate  models  which  work  well  in  a  priori  test 
but  which  fail  in  practice  due  to  the  inherent  instability  of  the  proposed  formulation.  The 
proposed  dissipation  model  has  been  incorporated  into  a  full  RST  closure  and  solved  for 
turbulent  channel  flow.  The  details  of  the  full  RST  closure  are  given  in  Natu39.  Results 
of  these  simulations  for  channel  flow  at  Re=590  are  show  in  Figure  13,  and  for  the 
rotating  channel  flow  case  in  figure  14.  These  full  model  results  are  highly  dependent  on 
the  chosen  pressure-strain  model.  They  are  not,  therefore,  an  indication  of  the  accuracy 
of  the  dissipation  tensor  model.  They  are,  however,  an  indication  of  the  dissipation 
tensor  model’s  stability  and  computability. 


6.  Conclusion 

The  proposed  model  for  the  dissipation  tensor  deviates  from  all  prior  models  in  its  use  of 
terms  which  involve  gradients  of  various  turbulence  quantities.  The  appearance  of  these 
gradient  terms  is  not  particularly  surprising.  There  is  no  reason  to  believe  that  the  source 
terms  in  the  RST  equations  should  not  be  functions  of  such  gradients  and,  in  fact,  there  is 
every  reason  to  believe  that  gradients  should  dominate  in  regions  of  strong 
inhomogeneity,  such  as  near  walls. 
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Figure  13.  Channel  flow  (Re=590).  Circles  denote  DNS  data  from  Ref.  36  and 
thick  line  denotes  full  RST  model  predictions. 


When  attempting  to  add  gradient  information  to  a  model,  the  variety  of  choices  is  so  large 
that  the  traditional  approach  of  expanding  a  quantity  (such  as  the  dissipation)  in  terms  of 
all  the  possible  unknowns  becomes  intractable.  This  paper  has  demonstrated  a  rational 
approach  to  deriving  the  gradient  terms  in  the  model.  The  resulting  terms  (Eqn  E.7)  do 


Figure  14.  Full  RST  prediction  of  the  mean  velocity  in  rotating  channel 
flow  (Ro=0.15,  Re=194).  Circles  denote  DNS  data  from  Ref.  38  and  thick 
line  denotes  the  model  prediction. 
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not  have  any  additional  model  constants,  but  do  tend  to  greatly  improve  the  model 
performance  in  regions  where  gradients  are  large. 

It  is  interesting  to  note  that  the  resulting  gradient  terms  do  have  associations  in  other 
contexts.  The  term  2v{KV2)  n{Kv2)  n  is  a  common  modification  of  the  dissipation  near  a 

wall  such  that  time  scales  at  the  wall  remain  finite,  and  the  term  vKn{^)  n  appears  in  the 
Reynolds  stress  anisotropy  transport  equation. 

The  quasi-homogeneous  term  in  Eqn  E.7  is  in  many  senses  the  harder  part  of  the 
dissipation  tensor  to  model.  Our  quasi-homogeneous  model  introduces  a  parameter  to 
represent  the  affects  of  mean  strain  on  the  dissipation.  In  theory,  this  term  should  only  be 
present  in  equilibrium  situations  (when  production  is  roughly  equal  to  dissipation).  The 
near-wall  2C  term  does  not  introduce  any  addition  model  constants  but  is,  like  the  strain 
dependent  term,  motivated  rather  than  derived.  Fortunately,  these  two  modifications  are 
not  large  in  most  flow  situations.  The  bulk  of  the  dissipation  model  is  carried  by  the  term 
g-fejSyK R.j .  This  model  satisfies  an  RDT  limit  and  is  realizable  in  the  2C  limit. 

It  agrees  reasonably  well  with  data  from  axisymmetric  expansion  at  both  small  and  large 
anisotropies.  The  key  innovation  in  this  term  of  the  model  is  the  functional  form  of  the 
blending  parameter/.  It  is  suggested  that  the  common  practice  of  expanding  parameters 
in  simple  polynomial  series  can  be  detrimental.  Such  expansions  do  not  perform  well 
when  the  expansion  variable  (such  as  F)  is  0(1). 
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Appendix  F:  Determination  of  Model  Constants 


To  callibrate  other  model  constants  like  ,  Cp2 ,  Cp2  and  C*2  ,we  used  steady  state  DNS 

data  for  channel  flow  (Re  =  590).  First  we  backed  out  the  constants  using  the  equation 
(14)  of  the  report  for  each  Reynold  stress  assuming  Cpl . 

Using  the  equation  (14)  for  R33 ,  we  get, 
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Using  equation  (F.l),  we  can  get  the  expected  value  of  Csp2.  After  that  if  we  put  this 
calculated  value  of  Cp2  in  equation  (14)  and  solve  the  equation  for  Ru  we  will  get  the 
expected  value  for  Cp2 . 
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Using  calculated  values  C*2andC”2  from  equation  (F.l)  and  (F.2)  and  then  solving  for 
Rn ,  we  get  expected  value  of  ,  given  by  the  equation  (F.3), 
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Then  we  tried  to  model  these  constants  using  different  scalar  quantities  like  F  and  a. 
First  we  matched  the  values  from  different  combination  of  these  scalar  varients  with  the 
expected  values  and  then  tryied  to  use  them  to  run  the  model.  After  doing  this  excercise 
for  number  combinations,  the  best  results  we  got  for  the  match  and  for  the  simulation  of 
channel  flow  (Re=590)  using  these  constant  in  the  model  are  shown  in  Figures  1,  2  and 
3. 

The  modeled  constants  are, 


Csp2  =-^ - 0.2F 

v+vT 

C'pl  =  0.2F2  -  0.006  £ 
£ 


Cw  = 

'~-p2 


Figure  1:  Backed  out  value  of  Csp2  from  DNS  data(channel  flow  590)  compared  with 
modeled  value 


Figure  2:  Backed  out  value  of  Cwpl  from  DNS  datafchannel  flow  590)  compared  with 
modeled  value 


Figure  3:  Backed  out  value  of  from  DNS  data(channel  flow  590)  compared  with 
modeled  value 
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Appendix  G:  PDF  ModdResultsmChaiindFtow 


The  simulation  of  the  proposed  RST  model  has  been  done  for  channel  flow  at  different 
Reynolds  numbers,  to  check  the  validity  of  the  proposed  model  for  dissipation  tensor  and 
pressure  strain.  The  following  figures  1  and  2  show  the  resulting  Reynold  stresses  from  the 
simulation  using  the  proposed  RST  model  till  we  get  a  steady  state  answer  for  the 
channel  flow  (Re  =  590) 


Figure  1  (a)Channel  Flow(Re  =  590)  Figure  l(b)Channel  Flow(Re  =  590) 


Figure  2  (a)Channel  Flow(Re  =  590)  Figure  2(b)Channel  Flow(Re  =  590) 


After  getting  R-tjl K,  we  simlply  multiplied  it  by  K  to  get  only  R{j  and  then  compared  it 
with  DNS  data.  Following  figures  3(a),  3(b),  4(a)  and  4(b)  show  the  prediction. 


Figure  3  (a)Channel  Flow(Re  =  590)  Figure  3  (b)Channel  Flow(Re  =  590) 


Figure  4  (a)Channel  Flow(Re  =  590)  Figure  4(b)Channel  Flow(Re  =  590) 


Figure  5  (a)Channel  Flow(Re  =  590)  Figure  5(b)Channel  Flow(Re  =  590) 
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As  can  be  seen  from  above  figures,  the  predictions  from  the  model  match  quite  accurately 
with  the  DNS  data.  All  the  Reynolds  stresses  are  modeled  very  well.  The  most  important 
prediction  for  a  turbulence  model  is  the  average  velocity.  Figure  6  shows  model 
predictions  of  horizontal  velocity  with  DNS  data. 


Figure  6  Channel  Flow(Re  =  590) 


We  also  tried  to  simulate  channel  flow  (Re=  395)  using  the  model,  to  check  the  model 
predictions.  Results  are  shown  below. 


Figure  7  (a)Channel  Flow(Re  =  395)  Figure  7(b)Channel  Flow(Re  =  395) 
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Figure  9  (a)Channel  Flow(Re  =  395)  Figure  9(b)Channel  Flow(Re  =  395) 


Figure  10  (a)Channel  Flow(Re  =  395)  Figure  10(b)Channel  Flow(Re  =  395) 
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Figure  1 1  (a)Channel  Flow(Re  =  395)  Figure  1  l(b)Channel  Flow(Re  =  395) 


Figure  12  Channel  Flow(Re  =  395) 
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Appendix  H:  PDF  ModdResuttsforH^ 


The  eddy  collision  model  has  been  applied  to  a  number  of  homogeneous  shear  flows. 
These  are  experimental  or  DNS  cases.  The  model  given  by  equation  1  predicts  all  but  the 
rotating  cases  reasonably  well. 


The  model  has  also  been  applied  to  turbulent  channel  flow  at  a  shear  stress  Reynolds 
number  of  590.  The  Reynolds  stresses  and  mean  flow  are  given  in  figure  5. 

The  fast  pressure-strain  constants  in  the  current  eddy  collision  model  are  given  by, 

Cpls  =  - 0.2F  Cplw  =  - 0.4F  CpV  =  0.2F2  -  0.006^ 

V  +  VT  V  +  VT  £ 

where  vt=1Cm-ej-  and  F  =  determinant Rt. Ik)  where  the  Reynolds  stress  tensor  is 
exactly  given  by  Rtj  =  jv-v'f  and  the  modified  dissipation  is  given  by 
£  =€- 2v  3(£'2)  3(£'2) .  The  transport  constant  is  CfI  =  0.15 . 


The  dissipation  constants  are  standard  values  with  modifications  for  low  Reynolds 
numbers  and  high  rotation  rates. 


CeI  =1.45  Cr  Ce2=Cr(1.83-0.166e  4w) 

1  +  0.21RO  Ro  =  i.(WW)2 
r  1  +  0.12RO  K  u  ,J 

Ce3=C,  (0.33  +  0.5f) 

In  the  oriented  model,  the  dissipation  transport  equation  is  not  necessary  (i 


-vk/z2). 
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Appendix  I:  BGK  moments 


Conservation  of  mass  (or  probability)  requires  that  the  integral  of  the  PDF  be  equivalent 
to  one  for  all  time.  This  means  that  the  integral  of  its  time  derivative  must  be  zero. 
Starting  from  the  BGK  model  for  the  time  derivative  of  distribution  gives  the  following 
expression. 


(i.i) 


Since  both  distributions  integrate  to  1,  we  can  see  that 

i|/3v  =  0  (1.2) 

Conservation  of  momentum  requires  that  no  mean  flow  be  created  by  the  relaxation 
process.  The  mean  velocity  of  the  flow  is  equivalent  to  the  first  moment  of  the  PDF.  By 
taking  the  integral  over  all  velocity  space,  we  can  show  that. 

(i-3) 

Using  the  fact  that  the  velocity  is  an  independent  variable  (from  time)  and  splitting  the 
velocity  into  its  mean  and  fluctuating  parts  gives. 

£  jv,/3v  =  Jv,./3v - Uj  j(jxKyV2e~^dv-  Jv’f(f^K)_3/2e“^avJ 


(1.4) 


By  definition  Jv,/3v  =  w, .  The  second  integral  on  the  right  hand  side  is  equal  to  1  and 

the  last  integral  is  zero  since  it  has  an  odd  integrand,  so  finally 

£u>  ~ui  -°}  =  °  (1.5) 

The  Reynolds  transport  equation  is  obtained  by  multiplying  the  PDF  relation  equation  by 
v'jV'j  and  then  integrating  over  all  velocity  space. 
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Since  |v'(.  v';.  fdv  =  Ry  by  definition,  and  the  last  integral  must  be  isotropic 
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Appendix  J:  Relaxation  Model  Moments: 

Here  we  verify  conservation  of  mass  for  the  relaxation  models  derived  in  Section  4  of 
Appendix  C.  The  method  is  the  same  as  before  starting  from  the  relaxation  model  for  the 
PDF. 

lf*=Jc„ (*(***)  ~e 

with  v =  v;  -m,  ,  v  \  =  V,.  -  Uj . 
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The  integrals  can  be  evaluated  to  give 
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Similarly,  to  verify  conservation  of  momentum  we  continue  as  follows. 
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b-/  Pv 


)V 


(J.5) 


C„¥rU-C 


3f 


du,- 

IT  ~  V'M  2K  **/  2AT  2^+(B 


J— t((k-«)2«,  +2  K-«X  +2ATu,  + jv»',/3v)  (J.6) 


Conservation  requires  the  above  equation  be  equal  to  zero,  this  implies  that 

(2 K  +  (u -m)2)(m,.  -u,)  =  2RJun  -  un)  +  fv)  v’„  v'„  /dv  (J.7) 

From  Appendix  I,  we  see  that  if/ is  Gaussian,  the  last  term  on  the  right  goes  to  zero,  and 
u.  =  u.  =  Uj  confirming  conservation  of  momentum  for  Gaussian  PDFs.  For  non- 
Gaussian  PDFs  the  above  relation  must  be  satisfied. 


The  Reynolds  stress  transport  equation  is  also  derived  similarly 


aw. 


fv'  v'  fc  v'  v’  \  l£-(±JlKY3l2e  *k  -  3e  v"v" _ f 

JV  iV  jlPV~  3,  -  J°MV /V  j{2K^37r^)  e  2K  2JC+(«-fi)2  7 


> 


(J.B) 


By  substituting  in  the  relations  v'f  =  v(.  -  m,  and  v\  =  v,  - m,  the  integrals  can  be  reduced. 

T=  C„*  J(f  >(*, -U,))(f,+  (SJ 


-c, 


3f 


A/  2iP  2/C+(m 


( («„  +  2v  («„-«„)  +  v>'„)  /3v 


(J.9) 


since  tcK)  312  e  4*"9v  =  0  (due  to  the  odd  integrand),  we  get 


1 


3fl‘n  fr'n 


-  c*  ~a)lR« + 2<“.  -«„)  Jv>'/ v'.  /3v+  Jv»'. v'.  /3v) 

The  first  integral  is  reduced  in  terms  of  ‘hats’. 

TT  =  Cu  #(f Hj  +(«,-  -«;•)) 

-  #2y+(U2  ((«  -  +  2(«„ -W„)jv V  ’ .  V  /dv  +  Jv  ’  V  V  V  fdv) 


(J.10) 


(J.ll) 


(J.12) 


(J.  13) 


To  ensure  the  correct  dissipation  of  energy,  we  require  that  the  model  satisfies  the 
equation^  =&  =  -€. 

-i  ^  =  «  =  -C„#(^+|(i!-«)2) 

+C„%— J— r  (  (k  -  S)2  JC  +  («„  -  «„  >  Jv  »•„  i  Jv  •, v v  v’„  /9V) 

This  can  be  simplified  to 
K+\{u-ii)2  +■ 5^-  = 

-ufK  +  (un  - un )  Jv v v fdv+\  Jv v V  v>'„  fdv) 

or  finally 

(1 +iH“  -«)2  +  35t)(2A: + <“  -  “>!)  = 

(„  _  sf + ii  Jv  ■  V  •,  V  fdv  +  Jf  Jv v ' »'.  v fiv 

For  an  elliptic  Gaussian  PDF  the  above  equation  reduces  to 

f +  W  =  7k  Jv  >  K  A  +  2*.  A  )  (J-15) 

or  alternatively 

f  =  +  ^  (J-16) 

This  can  be  rearranged  to  become 

ic„-M+3f  (I17) 


(J.14) 


So  if  we  assume  Gaussian  form  for  the  PDF,  the  Reynolds  transport  equation  can  be 
written  as  follows. 


W-i  +  m  (t'HHM  +  W  MR*.R»+2R,:R«) 

This  simplifies  to 

dRiJ  _  £  1  /  2  j>2  _  ft  RnR«j  \ 

3;  K  j_i+£aa5aa.  \  3  iV  uij  IXij  K  f 

K  2  K2 

Rearranging  into  the  classic  return  model  form  gives 


(J.19) 


2 


K  *mRm 


_ _ e_  n  e  K  2 k2  (  p  _1  \ _ 1 _ /?  /?  Jl - £ - I - R  R 

dt  “  K  y  /if  1_i.+5Da2aiB.VA0‘  3  £2  i-i+Mi  'Nwi^mn  3  i-i+Mi  ni  n 

K  2K2  K  2K2  K  IK2 

And  when  written  as  follows,  the  values  of  Cr  and  Cn  become  apparent. 


dfy  _ £_  n _ e_ 

dt  ~  K  *j  K 


( i+i—SnlSlL  V  .  I 

ill  z*i  i  4  -i  ( d  _ifcs.  )+-£ _ =i — /?  .ft .  -ft  .ft 

1  £  .Mm  +  3  i..a  +  *A  Ka>j  3  A uij )  T  K2  1-L+MmL  \  «'■  nj  m  m 

*  2*2  *  2Jf2  y  *  2*T2 


,  <J.20) 

I)  '(J-21) 
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Appendix  K:  Moments  of  a  Gaussian  PDF 

If  we  have  a  PDF,  f,  of  elliptic  Gaussian  shape  we  can  write  the  PDF  as 
/  =  [{27tf  det(/?nm)]"1/2e”^fi“v " v" 

Since  only  RtJ  is  a  function  of  space,  this  implies  that  the  derivative  satisfies 


(K.2) 


This  allows  us  to  write  the  third  moment  as  a  second  moment  and  then  apply  the  chain 
rule. 


Jv»'„  fdv  =  ■ -«*,  Jv V =  -S„  J(^- / (K.3) 

By  Gauss’s  divergence  theorem  the  first  integral  term  goes  to  zero,  since/  ->0  at  infinity. 
Differentiating  the  second  term  gives, 

Jv v'nv'mfdv  =  Rkm\f(v \  Snk+v\  Sik  pv  =  Rnm  Jv *,  fdv + Rim  Jv fdv  =  0  (K.4) 

which  is  what  we  would  expect  since  the  PDF  is  an  even  function  and  the  integrand  is  an 
odd  function  (cubic). 

The  expression  (K.2)  and  the  chain  rule  is  also  useful  for  evaluating  the  fourth  moment  of 
an  elliptic  Gaussian. 

J  v'mv»'i//v  =  -^.J  v v v ’  ^dv  =  Rkj  J / dv  (K.6) 

Taking  the  derivative  gives  the  fourth  moment  in  terms  of  the  second  moments 
=  ^  f  /  ( v 'm  v '  A  +  v  > )  A*  +  v 'm  v Snk )  dv  =  R^Ry  +  RniRmj  +  RmiRnJ  (K.7) 
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Appendix  L:  General  relaxation  moments: 

Here  we  verify  conservation  of  mass  for  the  more  general  relaxation  models  derived  in 
Section  6  of  Appendix  C.  The  method  is  the  same  as  before  starting  from  the  general 
relaxation  model  for  the  PDF,  Eq.  (C.33).  Conservation  of  mass  then  requires  that  the 
right  hand  side  of  the  zeroth  moment  be  equal  to  zero. 


-  CM  *  Jn,[(2*)3  d M-v** ‘ ■  av  -  C„  *  a- ' ) 

J* >'= c“  *] n«  “c"  J5'- v"» fSv 


(L.2) 


Substituting  with  v  \  =  v',+  n(.  -  ,  and  recalling  that  Jv'm  fdv  =  0  gives 
\Pv  =  CM  ^n,-Cw  J(v'm  V.+  K  ~K){un -Un))fdv 

J-j^V  =  CM  27  IT/,  -Cu  27  +  (Mm  “  Um  )(Un  ~Un))~Q 


0*3) 

0*4) 


The  momentum  equation  is  the  first  moment 

Tf-=  C„  *n,  |v,[(2»)3  det(4,)r,,a^1'-‘- 


du 


0*7) 


3v-c*'*f^Kt£wv'/3v  (L-5) 

r  =  -«.))K+(«,  -»,))v,/3v  (L.6) 

%  =  C„  #n,  i,  -  C„  J*  •-  V v  V  /3v + ) 

-C„  -^T— - - .n - rf-^mn  («„  —  W„  )  +  /?„„  f  M„  -  W„,  )  +  f  M„  —  M,„  )  (m„  -  )  n„  ) 

m  2K  l+(un -un )(«m -um )jfa~R^  \  mP\  n  n /  np  \  m  m  /  V  m  m)\  n  n)  p) 

Conservation  of  momentum  therefore  requires  that 
CM  #n„  M„  =  Cw  -5%-— - .n — rl  fv'  v'  v'_  fdv  +  Rm„un 

+  {Rmp  («»  ~l “n  )  +  Rnp  (“m  “ ' K  )  +  {Um  K  )  («»  ” ' «„  ) )} 

which  simplifies  to 

[i  +  («„  -  «„)(«„  -  =  TIT-  JV'm  V'„  v’p  fdv  +  UP 

+17^,;’  (Rn,p  («„  -«»)  +  ^K-«„))  +  T^m  («*  ~ )(«„“«„ )“P 

and 

(“,  "  UP  )[!  +  (“«  -  “n  )(«*  “  = 

T17 tf/m1  (  JV'm  V'n  V /^V  +  ^  ( Un  "««)  +  (»«  “  «« )) 


0*8) 


0*9) 


(L.10) 
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Appendix  M:  General  Reynolds  transport  Equation 

Below,  the  Reynolds  transport  equation  is  derived  for  the  general  relaxation  model. 


4  =  n„  Jv,  v',[(2^det(R„)]-'J^°'-"-3v-  3 


(M.l) 


Substitution  gives 


•  =  CM %nu  |(v’  v'j+  (u,  -ut)(uj  -  «.))[( 2tt)3 det (Rnm)TV2e 


"  2k >+(«„  -«•«  )K  -<  )ntti  J  1  J  m  n  ^ 


(M.2) 


+n«(*, -««,)(*, ■*.*'.  ■**)  (m.3) 

If  we  choose  un  =  un  and  use  the  general  form  of  the  Reynolds  stress  model, 
if  =  “ 2^(n/m  RmJ  +  UJm  Rmi ) ,  we  arrive  at  the  following  equation. 

-(nim Rmj + njm  Rmi)= cM n„ (4 + («,  - ul)(uj  - M .) jv-  v'.v’mv’„ /av)  (m.4) 

Which  gives  us  the  definition  of  R . 

4  =’A(^‘'nRmj+TljmRmi)-(ul-Ul)(uj-Uj)  +  ^-  Jv';  V'.  V*m  v'„  fdv  (M.5) 
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